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Abstract 


In this thesis we present the first numerical study of gravitational collapse in braneworld within 
the framework of the single brane model by Randall-Sundrum (RSII). We directly show that the 
evolutions of sufficiently strong initial data configurations result in black holes (BHs) with finite 
extension into the bulk. The extension changes from sphere to pancake (or cigar, as seen from a 
different perspective) as the size of BH increases. We find preliminary evidences that BHs of the 
same size generated from distinct initial data profiles are geometrically indistinguishable. As such, 
a no-hair theorem of BH (uniqueness of BH solution) is suggested to hold in the RSII spacetimes 
studied in this thesis—these spacetimes are axisymmetric without angular momentum and lion- 
gravitational charges. In particular, the BHs we obtained as the results of the dynamical system, 
are consistent with the ones previously obtained from a static vacuum system by Figueras and 
Wiseman. We also obtained some results in closed form without numerical computation such as 
the equality of ADM mass of the brane with the total mass of the braneworld. 

The calculation within the braneworld requires advances in the formalism of numerical relativity 
(NR). The regularity problem in previous numerical calculations in axisymmetric (and spherically 
symmetric) spacetimes, is actually associated with neither coordinate systems nor the machine pre¬ 
cision. The numerical calculation is regular in any coordinates, provided the fundamental variables 
(used in numerical calculations) are regular, and their asymptotic behaviours at the vicinity of the 
axis (or origin) are compatible with the finite difference scheme. The generalized harmonic (GH) 
formalism and the BSSN formalism for general relativity are developed to make them suitable for 
calculations in non-Cartesian coordinates under non-flat background. A conformal function of the 
metric is included into the GH formalism to simulate the braneworld. 



Preface 


Chapter 1 is the introduction and nothing in the chapter is original work. The author rederived 
the formulae in the literature while keeping d (the dimension of spacetime) and e (the sign to 
characterize spacelike or timelike foliation of spacetime) general. 

All the works in the rest chapters (Chap. 2, 3, 4 and 5), except for the works properly cited, 
are original. Within these: 

- A significant part of the work regarding initial data (presented in Chap. 4) was conducted 
in the collaboration with Evgeny Sorkin. Specifically, Evgeny proposed “direct method” 
(defined in the chapter), taught the author to calculate the total energy in general relativity 
following [87], and proposed to study the mass-area relation. 

- Eq. (5.34), the form of the relation among the circumstances of black holes, was proposed 
by Toby Wiseman. 

- All the remaining works were carried out by the author, with the guidance of Matthew W. 
Choptuik and William G. Unruh. 
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Please pay particular attention to sign conventions: the sign of metric, the sign of the Christoffel 
symbol, and the sign of the Riemann tensor. The conventions we employ in this thesis, are the 
same as those in Baumgarte-Shapiro [34] (therefore the same as those in Wald [32]). The sign of 
extrinsic curvature is the same as Baumgarte-Shapiro [34] (therefore the opposite of Wald [32]). 
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Chapter 1 


Introduction 

1.1 Overview 

Our observable universe is 3 +1 dimensional. The exploration on the possibility of the existence 
of extra dimensions can trace back to the work by Kaluza and Klein in 1920s [1, 2], where they 
tried to unify electromagnetism with gravity by using the metric tensor in five dimensional (5D) 
spacetime in which our universe was a 3+1 dimensional hypersurface. The basic idea of braneworld 
scenarios is that our observable universe could be a 3+1 dimensional hypersurface (the brane) 
embedded in a higher dimensional spacetime (the bulk). Only gravity can propagate freely into 
the bulk while all the matters are confined onto the brane. Among all the braneworld models, the 
single brane model constructed by Randall and Sundrum (also known as RSII) [3], is remarkable for 
its simplicity. It assumes one extra dimension with infinite size, a negative cosmological constant 
A in the bulk, tension A on brane which makes the brane a gravitating object, and Zi symmetry of 
the bulk with respect to the brane. A is set to a value such that general relativity (GR) is recovered 
on the brane at low energies [4]. However, the high energy behaviour, where the dynamics could 
be rather different from GR, is not understood as clearly. 

Black holes (BHs), among others, are objects that can illustrate the high energy behaviour. BHs 
are predicted by GR, and there is strong observational evidence for their existence [5], therefore 
RSII needs to reproduce BHs in order to be a viable physical theory. In the absence of matter, there 
exist a particular class of solutions in RSII constructed from vacuum solutions of four dimensional 
(4D) GR using the construction method in [6]. These solutions are called black strings when 
the corresponding 4D solutions are black holes. However, the black strings are unstable due 
to the Gregrory-Laflamme instability [7, 8]. In fact, Lehner-Pretorius [69] showed that certain 
type of black string (which is different from the black string in braneworld) evolves into a series 
of 3D spherical BHs connected by black strings and the black strings continue to evolve in the 
same pattern, which yields a self-similar configuration. Within finite asymptotic time, a naked 
singularity is created, which violates the cosmic censorship hypothesis. Besides, the Kretschmann 
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scalar R^vapR^ 1101 ^ (where R^ap is the Riemann tensor) diverges at the AdS horizon [6], which 
means the black string (in RSII) solution itself has pathology. Therefore, black strings can not be 
formed via natural processes such as gravitational collapse. 

What is the bulk counterpart of the 4D black hole? Noticing that the instability is most severe 
at the AdS horizon where black strings might “pinch off”, Chamblin-Hawking-Reall [6] proposed 
that BH with finite extension into the bulk was a stable state and was actually the end state of 
gravitational collapse. The investigation on this proposal started from studying one aspect of this 
proposal: the existence of such black objects. Tanaka [9] and Emparan et al. [10] conjectured, via 
holography, that static black holes with radii much greater than the AdS length, can not exist in 
RSII. Fitzpatrick et al. [11] and Gregory et al. [12] then argued that such a conjecture may not 
be justified since the arguments in [9, 10] did not take into account the strongly coupled nature of 
the holographic theory. 

Several other groups investigated this proposal using numerical techniques and tried to find 
specific solutions describing the BHs. As a first step, Shiromizu-Shibata [13] studied time sym¬ 
metric initial datasets for the Einstein’s equations, demanding that the data contained apparent 
horizons, so that subsequent evolution would presumably settle down to black holes. Kudoh- 
Tanaka-Nakamura [14] carried out calculations on static BHs in the braneworld, finding solutions 
with finite extension into the bulk, but only for radii that were small compared to the AdS ra¬ 
dius. Subsequent work by Kudoh [15] showed that corresponding solutions could also be found in 
the 6D braneworld, and that the horizon sizes could be larger than the AdS radius in that case. 
Further numerical calculations in 5D by Tanahashi-Tanaka [16] found time symmetric initial data 
with large apparent horizons. In 2008, Yoshino [17] repeated the calculation done in [14], using 
a more accurate numerical scheme, and argued that due to certain systematic errors exhibited by 
the solutions—and hence presumably affecting the previous computations—there was no evidence 
that BHs with finite extension into the bulk could exist at all (i.e. not even small ones). In 2011, 
Kleihaus et al [19] also claimed the same. 

Finally, the debate was settled by the remarkable work carried out by Figueras-Wiseman [20] in 
2011. Via perturbing AdS/CFT solution (which itself was also numerically constructed [21]), they 
obtained static black holes with a range of sizes, including ones large compared to the AdS scale. 
Recently (2014), Figueras presented an improved calculation [22] which obtained black holes with 
much wider range, thus presumably black holes can be of any size. In 2012, large BHs were also 
obtained by Abdolrahimi et al [24] using a different method, and there was preliminary evidence 
that their solution was the same as the Figueras-Wiseman solution. 
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Without proving the stability and uniqueness of the solution however, it is not clear whether the 
Figueras-Wiseman BHs can be formed via natural processes such as gravitational collapse. More 
generally , despite the success on the statics side, the full dynamics of RSII—a more important and 
interesting aspect—remains poorly understood, especially in the high energy regime. The answers 
to fundamental questions, such as the end state of gravitational collapse, or the interaction between 
the brane and bulk, remain unknown or vague. 

The main goal of this thesis work is to study the full dynamics in RSII via numerical calculation. 
Specifically, we will study the dynamical process of black hole formation as a result of gravitational 
collapse of massless scalar field (which lives only on the brane). To our knowledge, there is no 
previous calculations of the full dynamics of (any) braneworld scenario. The study is crucial to 
understand the dynamics of braneworld models. Technically, as it turns out, the study deepens 
our understanding and also extends the formalism of numerical relativity. 

The calculation in the braneworld is significantly more difficult than that in GR. The confine¬ 
ment of matter to the brane and the interaction between the brane and the bulk due to brane 
content including matter and brane tension, are the roots of the difficulties [25]. Gravitational col¬ 
lapse inevitably produces energies high enough to make the interaction between the brane and the 
bulk significant. This rules out any attempts to solve the problem using brane content only. Due 
to this interaction, the brane equations are not closed, thus the dynamics of the whole spacetime 
(including both the brane and the bulk) needs to be studied. 

Here we successfully performed a numerical study, and we directly show that the results of 
the gravitational collapse of sufficiently strong initial configurations are BHs with finite extension 
into the bulk. The extension changes from sphere to pancake (or cigar, seen from a different 
perspective) as the size increases. There are indications that BHs with the same size produced 
from different initial data, have the same shape, which means the details of initial data are lost in 
the resulting BHs. Therefore a no-hair theorem (uniqueness of BH) is suggested to hold also for 
BHs in the RSII spacetimes that are studied in this thesis. In particular, the BHs we obtained as 
the results of dynamical systems, are consistent with the ones previously obtained from a static 
problem by Figueras-Wiseman [20]. Please refer to Fig. 1.1 for a preview. 

There are gaps between the existing knowledge in numerical relativity and the knowledge needed 
for the numerical calculation of braneworlds. We developed some of these to make our numerical 
calculation possible. With these developments, we solved the regularity problem in numerical 
relativity (the irregularity associated with simulations in axisymmetric or spherically symmetric 
spacetime which appeared in previous numerical relativity calculations) by studying it from a 
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different perspective than previous works did. Our work deepens the understanding of the problem, 
and can potentially completely solve the problem. We generalized evolution formalisms (generalized 
harmonic and BSSN) of GR to make them suitable for simulations in non-Cartesian coordinate 
under non-flat background. A conformal function of the metric is included in the formalism to 
simulate the braneworld. The constraint violations found in the braneworld calculations, are cured 
by imposing the constraints at the brane. 



Figure 1.1: The area of bulk apparent horizon as a function of the areal radius of the horizon 
on the brane. The coupling strength of the brane tension is proportional to 1 /£ (where t is the 
AdS length), which is invisible to the BHs whose sizes (the areal radii on the brane) are much 
smaller than i. Therefore small BHs behave as 5D Schwarzschild. When the size is much larger 
than £, the asymptotic relation is that of the corresponding black string. Our data obtained from 
the evolutionary system, is consistent with the data obtained by Figueras-Wiseman [20] from a 
static system. 


1.1.1 Overview of Our Approach 

The system is formulated as an initial value problem (IVP), with Einstein’s equations in the 
bulk as the governing equation. The brane content, including brane tension and the matter, is 
imposed via Israel’s junction conditions [26], which serve as parts of the boundary conditions of 
the IVP. The matter on the brane evolves as ordinary 4D matter (which only “feels” the 4D metric 
on the brane). 

Given that calculations of this sort (or even published attempts) do not exist, and due to the 
prohibitive computational cost of performing the calculations in the fully 5D context, we restrict 
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ourselves to spherical symmetry on the brane, which makes the system axisymmetric in the bulk 
(i.e. functions depend on two spatial dimensions plus time). 

Even with this significant simplification, we are faced with a challenging task, not least since 
there are several features of the problem that have not been addressed, or fully resolved, in previous 
numerical relativity studies. Apriori, there are the following challenges 

(1) The numerical treatment of delta-function (distributional) matter 

(2) The numerical treatment of a brane with tension 

(3) Regularity issues induced by axisymmetry 

(4) Appropriate evolution schemes for use in non-Cartesian coordinate systems, and under non-fiat 
backgrounds 

(5) Properly incorporating the AdS boundary conditions 

(6) Higher dimensional numerical relativity 

Among the features mentioned above, some are more challenging and cause more severe prob¬ 
lems than others, and could be made into good projects by themselves. In fact, Oliver Rinne’s 
PhD thesis (at University of Cambridge, 2005) [71] was on regularity issues in axisymmetry. He 
studied a specific formalism (Z(2+l)+l system) for axisymmetric system. On the other hand, 
what we are going to present in this thesis is a “direction change” of the study on this topic. 
We will conjecture that the regularity problems are caused by the the fundamental variables used 
in numerical simulations, rather than coordinates. We further developed a few general methods 
to produce formalisms that yield regular results. Hans Bantilan’s PhD thesis (at Princeton Uni¬ 
versity, 2013) [68] was on the AdS spacetime in 5D, where lightlike signals can propagate to the 
spatial infinity and come back to its departure location within a finite proper time as measured 
by a timelike observer. In this thesis, by employing the formalisms we developed to simulate 
braneworld spacetime (in Chap. 3), there are no issues associated with the higher dimensions. The 
spacetime of the braneworld is asymptotically a part of the Poincare patch of the AdS spacetime 
whose causal structure is similar to that of Minkowski spacetime. Therefore Cauchy surfaces exist 
in the spacetime (so that Cauchy problem is well-defined), and the above mentioned problem as¬ 
sociated with lightlike signals will not occour. There will be complications related to the Killing 
horizons of the Poincare patch, however, such as a timelike curve of finite length can reach the 
Killing horizons. Please refer to Sec. 1.4 for more details about AdS spacetime, the Poincare patch 
and the background of the RSII braneworld. 
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During the investigation, we encountered peculiar properties and challenges which were not 
anticipated. In solving for time symmetric initial data, in most situations the problem can reduce 
to solving the Hamiltonian constraint for one unknown variable. In RSII, however, there have to be 
at least two unknown variables in the initial data to make Israel’s boundary conditions consistent; 
normally the shift and lapse functions can be arbitrarily specified utilizing the gauge freedom 
in choosing coordinates, while they have to satisfy certain conditions in the braneworld due to 
Israel’s boundary conditions. In most situations, the constraint violation problem in free evolution 
schemes, such as the generalized harmonic formalism and BSSN formalism, can be considered 
solved by constraint damping method [58, 61-63]. However, the same method does not control the 
severe violation at the brane of the braneworld. We found that we could solve the violation by 
properly imposing the constraints at the brane boundary. 

Therefore we add the following two additinal problems to the feature list 

(7) Peculiar properties in initial data problem 

(8) Constraint violation in braneworld 

In the following chapters, we will show how we solve these challenges and achieve the simulation 
of the dynamical process of black hole formation as the result of gravitational collapse on the brane. 
The thesis is divided into three parts. In the remaining sections of Chap. 1, we will introduce the 
braneworld in more detail plus some basic aspects of numerical computation. The second part 
is devoted to develop the machinery to perform simulations and to extract physical informations 
in the braneworld. This part consists of Chap. 2 and Chap. 3. Chap. 2 develops the conceptual 
aspects, such as the smoothness of apparent horizon across the brane, the relation between brane 
horizons and bulk horizons, the coordinate gauge at the brane, and energy aspects of braneworld. 
In Chap. 3, our approach to the regularity problem will be presented. We also show how evolution 
schemes of GR, such as the generalized harmonic formalism, can be developed and generalized in 
non-Cartesian coordinate under non-flat asymptotic spacetime background, to be made suitable 
for simulations in the braneworld. The final part is the numerical simulation. Initial data is 
obtained in Chap. 4, which also provides numerical results to support the discussion on energy 
aspects in Chap. 2. In Chap. 5, we show how the constraint violations are solved, and also show 
the simulation per se, and the physical results. Other than Chap. 1, all the works in this thesis, 
except for the ones properly cited, are original. 


6 




1.2. Spacetime Foliation and the Decomposition of Einstein’s Equations 


1.2 Spacetime Foliation and the Decomposition of 
Einstein’s Equations 

In this section, we will introduce the notations and conventions in the context of the (d— 1) + 1 
decomposition formalism of general relativity (GR), and Israel’s junction conditions, which are 
needed to introduce the braneworld. Here d is the dimension of the whole spacetime. Metric 
signature is (—, +), which is the same as that in [29, 32, 34, 36]. 

The (d — 1) +1 decomposition of Einstein’s equations, is a generalization of the 3+1 formalism 
of GR [34, 36]. We start with a (d — 1) dimensional hypersurface E embedded in d dimensional 
spacetime M. E can be thought as the brane in braneworld, or a “time”= constant hypersurface 
in an evolutionary problem. I.e. the formalism introduced in this section applies to both cases. 
Especially, the parameter t in this section, is not reserved for “time”. 

E divides M into two parts A'^. Please refer to Fig. 1.2. 



Figure 1.2: The embedding of a hypersurface in a higher dimensional spacetime. The whole 
spacetime M is composed of M + (the half sphere) and M~ (the plane with a big hole in the 
middle to fit M + ). The intersection of Al + and M~ is hypersurface E (the circle), ri^ are unit 
normal vectors of E, as seen in M ± . The convention is that n + points into M + , and n~ points 
out of Al~. Note: (i) the extrinsic curvature of E in M + can be different from that in M~. For 
example, in this figure E is a curved line in M _ , but a straight line in M + ; (ii) In general n + n~ 
(but a few authors incorrectly assumed n + = n~ when deriving junction conditions). For example, 
in the figure n + is not defined in M~ and n~ is not defined in A'I + . Another example is n~ = —n + 
in Randall-Sundrum braneworlds. 


Generally E can be a hypersurface within a hypersurface family that locally foliates the space- 
time. Let us use parameter t to characterize the foliation, where t = constant is a specific hy¬ 
persurface (Et) within the family. t(A'I + ) > t(M~) is the convention we adopt. Let n ± be unit 
normal vectors of E, as seen in M ± . The convention is that the direction of n ± is the same as the 
increasing direction of t. Therefore n + points into M + , and n~ points out of M~. The super-index 
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± is omitted whenever we express anything that is valid for both signs. n a is normalized as 

n a n a = e, (1.1) 

where e is +1 when E is the brane in the braneworld and —1 when S is a “time” =constant 
hypersurface. In terms of t, the normal vector n ;i is 

n,, = e oV,,/. (1.2) 

where a = (-^/eV^fV^f) 1 is the lapse function. V is the covariant derivative associated with the 
metric g The greek index p, v = 0,... (d — 1), where 0 is the time coordinate. Define 

= cm M , (1.3) 


which satisfies m^V^f = 1. 

The decomposition of Einstein’s equations is expressed in terms of the following quantities. 
The induced metric on E 

Ivv =9in' - e n»n v . (1.4) 

The extrinsic curvature is defined as 

= -l^lJ^ccnp = (1.5) 

Note the geometrical meaning: the extrinsic curvature is defined as the “change rate” of the normal 
vector n a along the hypersurface, as seen in the bulk (M ± ). Since the intrinsic observer on E can 
not even “feel” n a , the extrinsic curvature characterizes the extrinsic nature of the embedding. 1 
The definition of Riemann tensor is the same as that in [29, 32, 34, 36], which is defined via an 
arbitrary vector field v M 

- V^V a ^ = R\ a0 vr (1.6) 

We also adopt the same sign convention of Christoffel symbols as that in [29, 32, 34, 36] 

F /3 7 = M (fiV/3 ,7 “t" — 9Pi,a )) (1-^) 

r Thf‘ extrinsic curvature can also be obtained as the measure of the deviation of the geodesics in £ and A /, 
produced by the same vector lying within S. Please refer to section 2.3.3 and appendix B for the study from this 
perspective. 
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where = d M . 

1.2.1 The Decomposition of Einstein’s Equations 

Repeating the Gauss-Codazzi-Ricci decomposition [29, 32, 34, 36] while keeping d and e general, 
the Einstein tensor G= R MI/ — ^g^R is decomposed as 

{d) G^n v = i (-e • (d -% + K 2 - ; (1.8) 

( VY. = D « K - D l>K'a; (1.9) 

{d) G^r a Yp = ^ (A» K af} - la pC m K) + G-D Gaf) 

+ e (2K aii K^ - KK a p + )p a p ( K 2 + K^K^)) 

+ ^ (la/ 3 D^D^a- D a Dpa), (1.10) 

where D is the covariant derivative associated with the hypersurace metric 7 M „. £ m is the Lie 
derivative with respect to m M . The super-indices ( d ) and (d— 1) are here to characterize dimension. 
For example, ( d-1 )G Mv and < ' d ~ 1 (R are the intrinsic Einstein tensor and the intrinsic Ricci scalar of 
E. 

Imposing Einstein’s equations in the d dimensional spacetime (which relates the geometry with 
matter via matter’s stress tensor T MJ/ ) 


{d) G „„ = k d T ( 1 . 11 ) 

we obtain Hamiltonian constraint, momentum constraint and evolution equation 

k d p = | (-e^-VR + A' 2 - ; (1.12) 

ek d S a = D a K - D lt K^- (1.13) 

k d S a p = 1 (C m K aP - la pC m K) + (d -^G af) 

+ e - KI\ a p + \p a p ( K 2 + K^K^)\ 

+ i {{lapD^D^a - D a Dpot) , (1-14) 


where k d = 8i:G d and G d is Newton’s constant in d dimension. We have defined p = T^n^n 11 , S a = 
e T liu n u 'y ll a , S a p = T liv 'y ll a 'j v p, which yield the following decompostion: = pn^ + n^S v + 
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S^n,, + Sftv Taking trace gives T = S + ep, where S = S fiv ')^ v . 

The evolution equation (1.14) can be equivalently expressed as [34, 36] 


Evc\Hct(3 — cD^D 0 Q£ “t - OL 


-e {d - 1] R^ + KK a0 - 2 K a „K^ + e k d (S a0 - 


d -2 


(1.15) 


The definition equation of extrinsic curvature (1.5), can be equivalently expressed as 


£m7a/3 = —2aK a0 . (1.16) 

(1.15) or (1.14), together with (1.16) form a complete set of evolution equations [34, 36] using 7 M „ 
and as fundamental variables. This formalism is called ADM-York formalism of GR. 


1.2.2 Israel’s Junction Conditions 

In classical electromagnetism, there are situations where electric charges are highly concentrated 
on two dimensional surfaces such as the interfaces of different materials. The distribution of the 
electric charges is then singular in the 3D space. The electric fields appear to be discontinuous 
across the singular layers, and the relation relating the discontinuities of the electric fields with 
the singular distributions of the electric charges, is called a junction condition. 

Similarly, in GR, if the stress tensor is highly concentrated on E, it will cause discontinuities. 
The discontinuities are described by Israel’s junction conditions. 

Israel’s first junction condition [26] states that the intrinsic geometry of the hypersurface is well 
defined, i.e. the induced metric of E obtained from M + agrees with the induced metric obtained 
from M~ 

hnA = o, (i-i7) 

where the notation [a] = a + — a~. 

Integrating (1.14) over an infinitesimal layer across E, we get Israel’s second junction condition 

[26] 

kd or \ K uA = ek d ' l 1 - 18 ) 

where ,5^ I1V is the singular part of the projected energy-momentum tensor on E defined as 

,°+ 

.5% = / S lu Af (1.19) 

No¬ 
where d l = d t(dtY = adt is the proper length/time across E, and l is (arbitrarily) set to be 
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l = 0 at E. 

Integrating (1.12) and (1.13) gives 

rO + rO + 

y^ = / S M dZ = 0; o= I pdlmO, (1.20) 

Jo- J o- 

which is because the right hand sides (RHSs) of (1.12) and (1.13), as well as the non-£ m terms on 
the RHS of (1.14), are all well-defined and finite on both sides of the hypersurface. 

Eq. (1.13) and (1.18) give 

»„y" u = -Ki, (i.2i) 

which is the singular matter’s conservation law on E. 

Because of the identity [a 2 ] = [o]{a|, where {a} = a + + a - , eq. (1.12) gives 
kd\p] = \ ([A'j{A'} - j K^ v ){K^}Y Combining with (1.18), we have 

w= i yym - v<y = -yyn. a.^ 

Eq. (1.21) and (1-22) are constraint/conservation of the singular matter. 

1.2.3 Coordinate Description 

A coordinate system is needed to perform numerical relativity. A coordinate system using t as 
a coordinate, can be constructed as the following. On each E t , a coordinate system x l is assigned 
and x l is set to be differentiable across E t . Here the Latin index i runs from 1 to (d — 1). (t, x*) 
is then a coordinate system. 

Since (d t)^ V M t = 1 = m/*V M t, the shift function /3 M = {dtY — to m is perpendicular to V M t, 
therefore in the tangent space of E t . The evolution equations can be written in coordinate system 
by using £ m( * = C(g t y — Cp^, where £(a t )^ is simply dt in the coordinate system in which t serves 
as one of the coordinate [31]. In general /3 M = /?* (<9;) M and 7 M „ = 7 ^ (dx*)^ (dx-1)^. They can 
be further reduced to /?* and 7 ^ in this coordinate system by using ((?,:)'' = (0, 0,..., 1 ,..., 0) and 
(dx*) M = (0, 0,..., 1,..., 0), where only the i -th position is 1. An infinitesimal vector pointing from 
(t, x *) to (t + dt, x l + dx*) is then dt (dtY + dx* (<9,) M , and its length is derived as 

ds 2 = (e at 2 + PiP 1 ) dt 2 + 2/%dtdx* + 7 ydx*dx J ', (1-23) 

where Pi = jijP j . 
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1.3 Randall-Sundrum Braneworld II (RSII) 

In this section we will introduce the single brane model by Randall-Sundrum (RSII). We will 
first introduce the results regarding the general ideas of braneworld (there is one extra dimension, 
and matter is confined on the brane), then introduce the setup that is unique in RSII. 

From now on, both the brane and “time” =constant hypersurface will be discussed. Therefore 
we use different notations to distinguish between them. We will use n„, h^, 

( d_1 ^JZfjv to denote the extrinsic curvature, unit normal vector, induced metric, covariant deriva¬ 
tive, Einstein tensor, Ricci tensor of the brane (therefore = 1), and use K^„, n M , 7 MJ/ , 

D^, to denote the corresponding quantities of the “time”=const hypersurfaces 

(therefore n M n M = —1). Especially, the parameter t is reserved for “time”=const hypersurfaces. 

1.3.1 Formalism for General Braneworld 

A class of braneworlds are defined as a five dimensional spacetime where the matter is trapped 
on a (3+1) dimensional brane E, but gravity can propagate freely into the bulk. Since the brane 
E is our observed universe, we would expect that the induced equations on the brane become 
Einstein’s equations at a certain limit. Using Gauss’ equations and the d-dimensional Einstein’s 
equations, one gets [4] 

V-^Gap = k d |t '^h\h v p + h a p (t ^nV + | 

+ ()OC a p - /C Q7 /C%) - \h af} (1C 2 - K^K^) - E a0 , (1.24) 

where T a p is the stress tensor in d dimensional spacetime, and T is the trace. is defined as 

E^= {d) C a ^ 6 n a iEh\h s u . (1.25) 

is the Weyl tensor defined by the following relation [4] 

^ = l * + ( i _ ^ (9^ *■ ^Rs]ii ~ 9v[j * — Jd — l )(d. — 2) ^ ^ 9n['i9s]vi (1-26) 

where we used the notation A[ a( gj = ^(A a p — Ap a ). 

From (1.24) one sees that the matter content needs to be specified. The extrinsic curvature 
and the projection of the Weyl tensor need to be addressed. 
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1.3.2 Randall-Sundrum Braneworld II 

In RSII, the size of the extra dimension is infinite. The bulk is empty except for a negative 
cosmological constant. Therefore the stress tensor is 


T a p — —A g a p + y 2 / 38 ( 1 ), 


(1.27) 


where A is the cosmological constant in the bulk, l is the proper length discussed in section 1.2.2. 
The brane is located at l = 0. In general, l is only defined in a neighbourhood of the brane. The 
brane content is 

y a /3 = -Ah a/ 3 + r a p, (1-28) 


where T a p is the stress energy tensor of matter on brane, and A is the tension of the brane, which 
is required for the consistency of the theory (see below). 

There is a Z 2 symmetry with respect to the brane [3] in Randall-Sundrum braneworlds, so that 
m _ = — n + . Z 2 symmetry is mirror symmetry followed by an identification operation: a spacetime 
M with Z 2 symmetry can be obtained by taking a piece of spacetime with a boundary, then 
generating its image by parity transformation with respect to the boundary, and then gluing these 
two pieces together by identifying the piece with its image (i.e. identifying point p with its image 
as the same point, for all p G M). In another word, Randall-Sundrum braneworlds are “one-sided” 
spacetimes. This point will be made clearer in Sec. 1.3.4. As a consequence of n _ = —n + , one 
can obtain IC+g = — K,~g. Israel’s second junction condition (1.18) is then reduced to 


^ = -k d y aP - h. 




d -2 


= \ kd 


\ ha/3 , 

A —-2 + Ta/3 _ fta,S 


T 


d- 2 


(1.29) 


This equation relates JC a p with the matter distribution on the brane, which can eliminate the 
extrinsic curvatures in (1.24). The E a p term in (1.24) is related to the geometry of the bulk and 
can not be eliminated easily. However, as analysed in [4], this term is only important at high 
energy regime. 

Now, (1.24) can be evaluated on either side of the brane, which may be realized by performing 
the evaluation at l 7 ^ 0 and then taking l —> 0 , because it does not make sense to do calculation 
exactly on the brane due to the delta-functions. One bonus relation obtained from this procedure 
is [E a / 3 ] = 0, which is a consequence of the Z 2 symmetry. The Einstein tensor on the brane is 
therefore [4] 

(d_1) 0a/3 = -Ad-iha/3 + k d -lT a p + k 2 d TT a p - E±p, (1.30) 
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where 


A “-‘= + ,L311 
k *~' = kix 4WT Y (L32) 

TTa/3 = 2rT g^ _ h ^ T + l (hapT^T^ - 2T^Tg 7 ) . (1.33) 

This recovers Einstein’s equations on brane at low energy regime since E^ v is also negligible in 
the low energy regime [4]. However, the behaviour could be quite different from GR at high energy 
scheme. Furthermore, due to the interaction between the brane and the bulk, the equations on 
the brane do not form a close system. Therefore the attempts to solve the dynamics by evolving 
only the brane content, is ruled out. 


1.3.3 Parameter Setting 

In order to let kd-i have the correct sign, (1.32) requires that [4] 


A > 0. 


(1.34) 


According to (1.31), A^-i, the equivalent cosmological constant on the brane, can achieve any 
value by tuning the value of A. In the case A^-i is taken to be zero, (1.31) yields 


2 {d - 2) 


(1.35) 


where t is the AdS radius in the bulk, which relates to A as 


1 = 


(d — l)(rf — 2) 


2A 


(1.36) 


Eq. (1.32) and eq. (1.35) imply 


G d -1 


V(d ~ 3) 2 • 27rG d 

l 


Gd , 


(1.37) 


where G n = fc ra /87T is Newton’s gravitational constant in n dimension. This relation clearly shows 
that t is related to the relation between the Newton’s constant on the brane with that in the 
bulk. In the theory of RSII, (. is a freely adjustable parameter, whose value can be determined by 
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experiments. Taking d = 5, (1-37) reduces to G 4 = G 5 • \Z8irG$/£. Choosing the unit k§ = £ = 1, 
this equation implies /04 = 1 . 

1.3.4 Vacuum Solution 

There are a class of solutions when matter is absent [ 6 , 42-44] 

ds 2 = e~ 2 ^ e dx'Mx 1 ') + Ay 2 , (1.38) 

where y G (—00,00) is the extra-dimension and the brane is located at y = 0. x M is the coordinate 
in the 4D section. h ltl/ does not depend on y, and can be any 4D vacuum solution of GR [ 6 ]. 
Define z = £e v ^ when y > 0, and z = £e~ v ^ when y < 0. In coordinate space, this is a two-to-one 
mapping which maps ±y to z = texp(|j/|/t) G [£, 00). However, this “two-to-one” feature is only 
superficial, because the (x M ,y) and (x M ,— y) are actually the same physical points, according to 
the Z 2 symmetry. Under coordinates (x M ,,z), (1.38) is changed into 

£ 2 

ds 2 = — [(/i MJ/ dx M dx I/ ) + dz 2 ] , (1.39) 

where z > £ covers the whole physical spacetime, and the brane is located at z = £. 

The simplest case is to let h IJV be 4D Minkowski spacetime. The corresponding 5D space is a 
part of the Poincare patch of AdS (Anti-de Sitter) space, and can be regarded as the counter part 
of Minkowski spacetime in the braneworld (see Sec. 1.4). h Mt/ can also take black hole solutions, 
and the corresponding solutions in the braneworld are called black strings. 

1.3.5 Matter in RSII 

Because of the specific form of the RS braneworld II, equation (1.21) and (1-22) now reduce to 

= 0; (1.40) 

0=0, (1.41) 

where D is the covariant derivatives associated with h Eq. (1.40) is important since it is the 

conservation law of the matter on brane. Since the tension part Asatisfies the conservation 
law, then eq. (1-40) requires that the matter part T a p must satisfy the conservation law. This is 
consistent with equation of motion of matter on brane, which takes its form in 4D GR since matter 
is strictly trapped on the brane and can not directly “feel” the extra dimension. 
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We will study the gravitational collapse of massless scalar field on the brane. For massless 
scalar field (denoted by $), the conservation law is equivalent to its equation of motion 

= 0. (1.42) 

The matter’s energy momentum tensor is 

Tftv = (1.43) 

1.4 The Global Structures of AdS Spacetime and the 
Poincare Patch 

As discussed in Sec. 1.3.4, in the context of the braneworld, the counter part of the Minkowski 
spacetime is 

ds 2 = — ^ — df 2 + dr 2 + r 2 (d# 2 + sin 2 9d(j) 2 ) + dz 2 J , where z > i. (1-44) 

The spacetime with this metric and z > 0, is a Poincare patch [75-77] of the AdS spacetime with 
AdS radius i. In this section we will discuss the global structure of AdS spacetime and its Poincare 
patch. The discussion is necessary for us to examine whether the RSII braneworld can be defined 
as an initial value problem, and whether event horizon can be defined. These two aspects are going 
to be discussed in Sec. 2.3.1. 

1.4.1 AdS Spacetime and the Poincare Patch 

The d dimensional homogeneous isotropic spacetime satisfying Einstein’s equations with a 
positive (or negative) cosmological constant, is called a de Sitter (or Anti-de Sitter) spacetime, 
and is denoted as dSd (or AdSd). In this section, we focus on the Anti-de Sitter spacetime. As 
discussed in [75-77], AdS^ is a hyperboloid of radius £ 

-X 2 -X 2 + j2Y 2 = -£ 2 , (1.45) 

i= 1 
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embedded in the fiat spacetime expanded by (Xq. Xi, Y\. whose metric is 


d—l 

ds 2 = -dA 2 - dX'l + ^ d^ 2 . (1.46) 

i= 1 

One coordinate system of the AdS spacetime is defined as 

Xq = l secy cos T, (1-47) 

Xi = £ secysinT, (1-48) 

l)=.£uytany, for i — 1, d — 1. (1-49) 


Here ay satisfy 


< 2-1 


I >* 2 = L 

i= 1 

i.e. they form a (d— 2) dimensional unit sphere in ( d— 1) dimensional Euclidean space. Under this 
coordinate system, the AdS spacetime is 


ds 2 


— (-dT 2 + dy 2 + sin 2 ydf^_ 2 ) , 


(1.50) 


where dfl 2 _ 2 stands for the line element on a ( d— 2) dimensional unit sphere in ( d— 1) dimensional 
Euclidean space. Eq. (1.50) with y £ [0, 7 t/ 2) and T £ (—oo,oo) defines the whole AdS spacetime, 
and is therefore a global cover of the AdS spacetime [75-77]. 2 
Alternatively, another coordinate system is 


*=u, 

Y i+1 = —, for i = 1,..., d — 2. 

z 


(1.51) 

(1.52) 

(1.53) 

(1.54) 


The range of the variables are z £ (0, oo), t £ (—oo, oo). X\, X2,Xd-2 form a (d— 2) dimensional 

2 In eq. (1-47) and (1.48), T is periodic with the period 2ir , therefore there exist closed timelike curves. The 
global AdS spacetime is to unfold these closed curves by removing the identification of T with T + 2tt. 
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Euclidean space. The metric of the AdS spacetime is now 


d—2 


ds 2 = — -dt 2 + dz 2 + ^ dx 2 


(1.55) 


i= 1 


Expressing the space expanded by Xj’s under spherical coordinates, this metric is 


ds 2 = — (-df 2 + dz 2 + dr 2 + r 2 d0 2 _ 3 ) 


(1.56) 


where dfl^_ 3 stands for the line element on a (d — 3) dimensional unit sphere in (d — 2) dimen¬ 
sional Euclidean space. This spacetime is called a Poincare patch of the AdS spacetime, and the 
coordinates (1.55) or (1.56) are called Poincare coordinates. 3 


1.4.2 The Penrose Diagram of the AdS 2 Spacetime 

In this section we mainly focus on the causal structure of the spacetimes, which can be con¬ 
veniently studied via the tool known as the Penrose diagram. For simplicity, let us start with 
AdS 2 —the 2D AdS spacetime. The metric of the global cover is 

ds 2 = (—dT 2 + d X 2 ) , (1.57) 

cos^ x 

and the ranges of the two variables are T £ (—oo, oc) and y £ (—7r/2, 7 t/ 2). 4 The Penrose diagram 
is shown in Fig. 1.3(a), where the T = const and y = const curves are simply horizontal and vertical 
lines. Eq. (1.57) shows that the null curves are straight lines with the slope dT/dy = ±1 (e.g. the 
orange lines in the hgure). 

An interesting feature of the AdS spacetime is that the lightlike signals can propagate to spatial 
infinities (y = ±7 t/ 2), and then come back to its departure location within a hnite proper time 
(measured at the departure location). Or precisely speaking, there exist closed causal curves with 
hnite proper length, which connect “local” regions with spatial infinities. As an example, let us 
consider a lightlike signal departing at point A in Fig. 1.3(a) and propagating along the orange 
line AB to point B (spatial infinity). It then propagates back from B to C along the orange line 
BC. For an observer sitting at y = 0, its proper time lapses 7r, which is hnite. On the other hand, 

3 (1.56) remains unchanged by the scaling: (t, z,r) —t ( ct,cz,cr ) where c is an arbitrary positive constant, and 
it can be either dimensionful or dimensionless. In particular, when taking c = 1/i , the (t, z,r) can be regarded as 
parameters with length dimension measured by unit £. In this section (Sec. 1.4) these coordinates will be treated 
as dimensionless parameters since it is easier to relate to the global cover. In all the other parts of the thesis, we 
will treat them as parameters with length dimension measured by unit t. 

4 Or equivalently but more closely to relate to (1.50) and higher dimensional AdS spacetime, we can also take 
X £ [0,7t/ 2) with "0 — 0,7r, where -0 is to parametrize in (1.49) as = cos0. 
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the proper distance from the observer to x = 7r /2 in a T = 


t/2 


1 


/ 0 cosy 


d% = 2 arctanh tan 


t/2 


const hypersurface is 
= 2 arctanh(l) = oo. 


(1.58) 


Actually the lengths of all the spacelike curves connecting a point at x = tt/ 2 with a point at 
X 7 t/ 2 are infinite, i.e. x = 7r /2 is truly a spatial infinity. For higher dimensional AdS spacetimes 
with d > 2, the above analysis is also valid. 

The metric of the Poincare patch under the Poincare coordinates is now 

ds 2 = — (—dt 2 + d£ 2 ) , (1.59) 


where t G (—oo,+oo) and z G (0,+oo). To get the Penrose diagram of the Poincare patch, we 
relate the coordinates (T, x, •■■) with the coordinates (t, z, r, ...) by the equality of eq. (1.47~1.49) 
and eq. (1.51~1.54). The coordinates (T, x, •••) can then be expressed in terms of (t, z, r,...). i.e. for 
every (t, z, r, ...), there is a corresponding (T, x, ■■•), which is a point on the Penrose diagram of the 
global AdS spacetime. The Penrose diagram of the Poincare patch, is then a subset of the Penrose 
diagram of the global AdS spacetime. The Penrose diagram for the Poincare patch is obtained as 
Fig. 1.3(b,c). In the figures, the blue lines are t = const, and the black lines are z = const. For 
the Penrose diagram, the boundary is divided into 

(1) Point i° represents z = oo, the spatial infinity. 

(2) Point i + represents t = +oc, the future timelike infinity. 

(3) Point i~ represents t = —oo, the past timelike infinity. 

(4) The line jf represents the future Poincare horizon (see, e.g. [39]), which is a Killing horizon 
associated with the Killing vector dt- As explained below in Sec. 1.4.3, this boundary can be 
reached by timelike curves of finite proper length, therefore these are not infinities. However, if 
the future oriented null curves are parametrized by the Killing parameter t, the null curves will 
end at this boundary when taking t —> oo. For example, the curve defined by (t, z) = (t,t+ 1) 
is a null curve. When t +oo, this curve ends at jf. i.e. this boundary appears to be infinity 
as measured by the Killing parameter t. In this sense, this boundary is an analogy of the null 
infinity of Minkowski spacetime [32]. We will put the word “infinity” into double quotes if it 
is not a true infinity but appears to be infinity as measured by t. 

(5) Similarly, j~ is the past Poincare horizon. 


19 




1.4. The Global Structures of AdS Spacetime and the Poincare Patch 


(6) z = 0 (which is also \ = ~ 7r / 2), the conformal boundary (see, e.g. [40]) of the Poincare patch, 
which is a part of the timelike spatial infinity of the global cover. 





Figure 1.3: Penrose diagrams of AdS 2 and its Poincare patch, (a) is the Penrose diagram of the 
whole AdS 2 - The T = const lines are horizontal, and \ = const lines are vertical. The range of 
the two variables are T G (—oo, oo) and \ G (—7t/2, 7t/ 2). Lightlike curves are locally straight lines 
with slope dT/d\ = ±1. (b) is the Penrose diagram of the Poincare patch, while (c) emphasizes 
that it is embedded into the global AdS spacetime. 


Fig. 1.3(b,c) show that in the Poincare patch there still exist lightlike signals going to spatial 
infinity (at the conformal boundary) and coming back whose trajactory can be connected by a 
timelike curve with finite proper time. e.g. the trajectory C —> A —> B in Fig. 1.3(c) is connected 
by the purple line BC whose proper time is finite. 

1.4.3 The Infinities of the Poincare Patch 

In the last subsection, we show that at the i 0 ,^ and jf, at least one of t = oo and z = oo is 
taken. Now we examine whether these are 
proper length of the t = 0 curve is 

d£ 

z 



“true” infinities, by evaluating the proper lengths. The 


f 1 dz f°° dz 

Jo z + Ji z 1 


(1.60) 
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where both terms on the right hand side are +oc. Actually the lengths of all the spacelike curves 
connecting a point at 2 = 00 with a point at z 7 ^ 00 are infinite. Thus i° is a true spatial infinity. 
Similarly, are true timelike infinities. 

For the Poincare horizons which appear to be null “infinities” measured in terms of the Killing 
parameter t. it is easier to perform the calculation under (T, y) coordinates. One can see that 
the proper time of FG —a curve connecting the past Poincare horizon with the future Poincare 
horizon as shown in Fig. 1.3(c)—is finite, which also means the lightlike signal propagating along 
F —* O —> G enters from the past Poincare horizon and escapes via the future Poincare horizon, 
within a finite proper time measured by an observer (whose world-line is FG). This means the 
information or entities in the space can “vanish” within a finite time. Also, the proper length of 
spacelike curve ED is finite. Therefore jf- are not true infinities. 

On the other hand, however, the Poincare patch (1.55) has a timelike Killing vector dt, and 
jf (the Poincare horizons [39]) are the Killing horizons associated with this Killing vector. In the 
literature, z = 00 is often referred as AdS horizon [ 6 ]. When approaching jf , the time parameter 
associated with the Killing vector (which is t) approaches infinity. An observer whose world-line 
is the integral curve of the Killing vector (i.e. his spatial coordinates of the Poincare coordinates 
remain constants), is called a Killing observer. For a Killing observer, his proper time takes 
(—(X), 00 ). is the future boundary of the spacetime region that can possibly be observed by any 
Killing observer, and j~ is the past boundary of the spacetime region that can possibly be affected 
by any Killing observer. 

1.4.4 The Poincare Patch of AdS$ 

AdS 2 catches the most important features of the more general AdS^- Also, AdS^ has a transla¬ 
tional symmetry in the space expanded by Xi s. Therefore, any finite x^s can be brought to Xi = 0 
using this symmetry, where the above discussions regarding AdS 2 still apply. On the other hand, 
however, the behaviour by “turning on” the space expanded by x^ s, might introduce some new 
features, because the Xi can take infinities. Therefore AdSs is a better model for the more general 
AdSd with d > 3. 

For AdS 3 , the metric of the global cover is 

£2 

ds 2 =- 5 — (—dT 2 -(- dy 2 + sin 2 yd?/> 2 ) , (1.61) 

cos^ y 
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where x € [0, 7 t/2 ), ip G [0, 27r) and T G (—oo, oo). The metric of the Poincare patch is 


ds 2 = — (-dt 2 + dz 2 + dr 2 ) 


(1.62) 


where z G (0, oo), t G (— 00 , 00 ) and r G (— 00 , 00 ). 5 Eq. (1.47-1.49) and eq. (1.51-1.54) now reads 


t ( d ~ 2 . \ 

X 0 = i secycos T = — 1 y~)x 2 — t 2 + z 2 + 1 I , 


11 

Xi = t secxsinT = —, 

z 


(1.63) 

(1.64) 


Y± = i u>i tan x — l cos ip tan x 
I 2 = £ 0 J 2 tan x = £ sin ip tan x 


£ r 


(1.65) 

( 1 . 66 ) 


where we have parametrized w»’s by io\ = cos ip and u >2 = sin ip. 

The Penrose diagram of the global cover is a cylinder, where the radius of the cylinder is x, the 
angle is ip and the vertical axis is T . Similar to AdS 2 , the Penrose diagram of the Poincare patch 
can be embedded into that of the global cover. Expressing ( T,x,ip ) by ( t,z,r ) via eq. (1.63-1.66), 
and then plotting in terms of (i, z, r), gives the Penrose diagram of the Poincare patch embedded in 
the Penrose diagram of the global cover. Please refer to Fig. 1.4. Fig. 1.4(a) is the Penrose diagram 
of AdS 3 and the Penrose diagram of its Poincare patch with some special lines. The vertical black 
line is the line of (x, ip) = ( 7 r/ 2 , it) which represents ( t , z, r) = ( t , 0, 0) with t G (— 00 , 00 ). The other 
two straight lines (which are half blue half black), correspond to the null “infinities” for r = const 
hypersurfaces. These three lines and the plane surrounded by them, form the Penrose digram of 
the Poincare patch of the AdS 2 discussed in Sec. 1.4.2. The three vertices of this triangle are i° 
and which represent the spatial infinity and timelike infinities, respectively. The boundary of 
AdS 3 is divided into (see Fig. 1.4) 

(1) z°, point (T, Xj ip) = (0, 7 t/ 2 , 0), the spatial infinity of z = 00 and r = ± 00 . 

(2) i + , point (' T,x,ip) = (zr, 7 t/ 2 , n), the future timelike infinity. 

(3) i~ , point (T, x,V0 = (— 7 t, 7 t/ 2 , 7 t), the past timelike infinity. 

(4) straight line described by (ip = 0 or 7 r, T = — xcosip + ir/2), the future null “infinity” with 
finite r. 

5 Again, equivalently but more closely related to higher dimensional spacetime, we can also take r G [0, 00 ) with 
6 = 0 , 7 r where 013 is parametrized as cos = cos 0. 


22 




1.4. The Global Structures of AdS Spacetime and the Poincare Patch 


(5) j z , straight line described by (ip = 0 or n, T = xcosip — 7 t/2 ), the past null “infinity” with 
Hnite r. 

The subscript z is to label the fact that the null “infinities” are within the t — z plane (since 
r = const). Similarly, the two crossing lines (half red half black) on the surface of the cylinder, are 
the null infinities for z = const surfaces, and are denoted as j ±, which are true infinities, i.e. 

( 6 ) j+, straight line on the side surface of the cylinder described by (y = 7 t/ 2 , T = ±ip with 
T > 0), the future null infinity with hnite z. 

(7) j~, straight line on the side surface of the cylinder described by (y = 7 r/ 2 , T = ±i/> with 
T < 0), the past null infinity with hnite z. 

( 8 ) z = 0 , the wedge portion of the side surface of the cylinder, the conformal boundary of the 
Poincare patch, which is a part of the spatial infinity (y = 7 t/ 2 ) of the global cover. 

The side surface of the cylinder corresonds to y = n/2. Its gray portion (on Fig. 1.4(a)) is the 
z = 0 hypersurface in the Poincare patch. Unfolding the cylinder along ip = 0, we get Fig. 1.4(b). 
The figure shows that the two crossing lines on the surface of the cylinder are actually straight 
lines with slope dT /dip = ±1, which are null curves: taking y = 7 r/ 2 , eq. (1.61) implies straight 
lines with dT/dip = ±1 are null. Fig. 1.4(c) emphasizes the null “infinities” of the Poincare patch. 
We take null curves parametrized by the Killing parameter t, and these null curves end at the null 
“infinities” when t —> oc. The orange lines are null curves with hnite r, which end at jf. The 
green lines are null curves with hnite z, which end at (and start from j~). At (or j±), z (or 
r) is hnite. On the other hand, the purple lines are null curves which end at the surface where all 
of f, r and z are inhnite. Therefore these null “infinities” with all of t, r and z approaching their 
infinities are denoted as 

(9) j+., the future null “infinity” with inhnite r and inhnite z. 

( 10 ) j~ z , the past null “infinity” with inhnite r and inhnite z. 

Within these null “infinities”, are not true inhnities, as discussed above for AdS 2 . Similarly, 

j-P z are not true inhnities either. On the other hand, i°, and j^r are located at the surface of 
the cylinder (where y = 7 r/ 2 ), and are true inhnities. The boundaries that can be reached by null 
curves are defined as the future/past Poincare horizons [39] 

( L67 ) 
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are the upper/lower boundary surfaces of the Penrose diagram (of the Poincare patch). The 
expression for these two surfaces can be obtained as follows [78]. Using eq. (1.63-1.66), ( t,z,r) can 
be expressed in terms of (T, x, ip) as [78] 


t sin T 

cos T — sin x cos ip ’ 

( 1 . 68 ) 

cosy 

z _ 

cos T — sin x cos ip ’ 

(1.69) 

sin x sin ip 

r =- 

cos T — sin x cos ip 

(1.70) 


Because z = oo at j , the equation for surfaces is described by requiring the denominator of 
(1.69) to be zero, which is 

cosT — sinycosV' = 0. (1.71) 

The upper/lower surfaces of Fig. 1.4(a,c) are generated using this equation. The upper surface is 
confirmed to be the future null “infinity” by Fig. 1.4(c) where the future null curves end at the 
upper surface. By the time reversal symmetry of the spacetime, the lower surface is also confirmed 
to be the past null “infinity”. 

Now we prove that the Poincare horizons described by (1.71) are indeed the Killing horizons 
for Killing vector dt- We define 

/ = cosT — sin y cos V 1 , (1-72) 

then the Poincare horizons are described by / = 0. For these hypersurfaces to be the Killing 
horizons, we need to prove [38] (1) the Killing vectors dt are orthogonal to these hypersurfaces; 
(2) the Killing vectors are null on these hypersurfaces. 

In general, the hypersurfaces described by / = const is orthogonal to the vector field d a f. d a f 
is evaluated as 

d a f = (— sinT, — cosy cos-0, sin x sin ip), (1.73) 

under the coordinates (T,x,ip)- A direct evaluation shows that d a f become null at the hypersur¬ 
faces specified by f = 0. On the other hand, a long but otherwise direct calculation shows that at 
the hypersurfaces specified by f = 0 

{d t ) a = (sin T/ cos 2 x) • d a f. (1.74) 

i.e. at the Poincare horizons, dt is proportional to d a f. Therefore the two conditions for the 
Poincare horizons being Killing horizons are met. 
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Figure 1.4: Penrose diagrams of AdS 3 and its Poincare patch. In (a), the cylinder is the Penrose 
diagram of the whole AdS 3 , while the portion surrounded by the gray surfaces (with different 
darkness) is the Penrose diagram for the Poincare patch. For the notations on the diagram: i° and 
are for points; j^r and jf 1 are for lines; j± z are for surfaces; and CB stands for the conformal 
boundary which is the wedge portion of the side surface of the cylinder. Note: j± z appearing on the 
side of the vertical black line, are actually for the surfaces rather than for the vertical black line, 
(b) is the conformal boundary (the z = 0 hypersurface of the Poincare patch), which is the part of 
the Penrose diagram on the side surface of the cylinder, (c) is to emphasize the null “infinities”, 
where the orange lines are the null curves with finite r, the green lines are the null curves with 
finite z and the purple lines are the null curves with all of t, r, z approaching their infinities. When 
taking t —> +oo, the orange lines end at j +, the green lines end at y+ and the purple lines end at 

7+ 

Jrz 
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The upper/lower boundary surfaces have another interesting property. We can put a second 
patch on top of the existing patch, by requiring that the i° of the new patch is the i + of the 
old patch. The new patch can be regarded as a transformation from the old patch, where the 
transformation is composed of the following two steps: (1) the patch is rotated by 7T with \ = 0 as 
the rotational axis. This operation is ip —> ip + n. (2) the patch is moved up vertically by 7r, which 
is T — > T + n. After these two operations, eq. (1.71) remains unchanged, i.e. it also describes the 
upper/lower boundaries of the new patch. Therefore, there is no gap between the two patches. 

1.4.5 The RSII Braneworld and Its Structure 

To get a better idea about the structure of the Poincare patch and the brane, we show the 
following hypersurfaces with one of the coordinates being constants. Fig. 1.5 shows the r = const 
hypersurfaces, where the blue lines are t = const. Or in another word, these blue lines are the 
integral curves of d z , and the black lines are that of <9*. Fig. 1.6 shows the t = const hypersurfaces, 
where the blue lines are the integral curves of d z , and the red lines are that of d r . Fig. 1.7 shows 
the z = const hypersurfaces, where the red lines are the integral curves of d r , and the black lines 
are that of <9 t . In particular, the brane in the RSII braneworld is the z = 1 hypersurface, which is 
shown in Fig. 1.7(b). 

The following discussion applies to general d case. 

As introduced in Sec.1.3.4, the spacetime background of RSII braneworld is to take the z > 1 
portion of the Poincare patch (which is called the bulk). Therefore the z = 0 boundary is eliminated 
from the Penrose diagram, and the global causal structure is the same as that of Minkowski 
spacetime. 

In fact, any z = zq (where the constant zo G (0,oo)) can serve as the brane, because the 
extrinsic curvature of the z = zo hypersurface can be calculated as 1 C MI/ = h^, which satisfies 
the Israel’s junction condition (1.29) applied to the vacuum case. Here h is the reduced metric 
on the hypersurface, and the extrinsic curvature is calculated based on the unit normal vector of 
the hypersurface, i/, pointing into the bulk (larger z direction). Within the AdS spacetime, the 
z = zo hypersurfaces are not geodesic surfaces (in the sense that the extrinsic curvatures of the 
hypersurfaces do not vanish, see Sec. 2.3.3), but are surfaces with constant acceleration, in the 
sense that every freely moving massive particle within & z = Zo surface has a constant acceleration 
as measured by the co-moving inertial observer in the AdS spacetime. Let d be the dimension of 
the AdS spacetime, and let be the d- velocity of a massive particle freely moving within a z = zo 
hypersurface. The path of the particle is then described by the timelike geodesics generated by the 
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tangent vector based on the intrinsic metric of the z = Zq hypersurface, which is v a T> a i= 0, 
where T> is the covariant derivative associated with h fiv . Let be the d -velocity of an observer, 
then the acceleration of the particle observed by this observer is 

= zt“V Q K - u»). (1.75) 


If the observer is an inertial observer in the AdS spacetime, his trajectory is then described by 
a timelike geodesics as rt“V a w A1 = 0. Furthermore, if the observer’s instant velocity is he is 
a co-moving inertial observer. The acceleration of the massive particle moving along a timelike 
geodesics in the z = zo hypersurface, as measured by a co-moving observer in the AdS spacetime, 
is then v “Vc/, which is described by eq. (2.27) (see Sec. 2.3.4) 

v a W a i = n^v a i>^lC a p = m >1 v a v^h a p = —n M . (1.76) 


This acceleration is a constant vector (as measured by the comoving inertial observer in the AdS 
spacetime) pointing out of the bulk. 



Figure 1.5: r = const hypersurfaces, where (a) is r = 0, which is the AdS 2 diagram, (b) is r = — 1 
and (c) is r = —8. The blue lines are the integral curves of d z , and the black lines are that of dt- 
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Figure 1.6: t = const hypersurfaces, where (a) is t ~ 0. (b) is t = 2 and (c) is t = 8. The blue 
lines are the integral curves of d z , and the red lines are that of d r . 



Figure 1.7: z = const hypersurfaces, where (a) is z = 0, which is opened as Fig. 1.4(b). (b) is 
z = 1, which is the brane in the RSII braneworld. (c) is z = 7. The red lines are the integral 
curves of d T , and the black lines are that of dt- 


28 






1.5. Evolution Schemes 


1.5 Evolution Schemes 

There exist various formalisms of GR, among which only the ones that are strongly hyperbolic 
(refer to, for example [34]), can be used as a well-defined formalism of an initial value problem. 
In this section we only (very) briefly sketch the generalized harmonic formalism since further 
developments will be present in the next chapter with more details. 

The generalized harmonic (GH) formalism [61] uses the gauge source functions 

H a ~ V^V/jx" = = -T a , (1.77) 

as fundamental variables. The notation ~ means the equation is a constraint relation. Einstein’s 
equations can now be written as 

1 9/j.iy,a0 ~ 9 P( t/ j,9v)p,a ~ ~ ^ \ta = ^ _ ^g^ v T\ . (1.78) 

A coordinate gauge choice can now be realized via specifying the i? M ’s. As long as H M does not 
include derivative of metric functions, the principal part of the above equation — ^g a ^g^,ag is 
manifestly strongly hyperbolic. 

Both the generalized harmonic formalism and the BSSN formalism [34, 45, 46] are widely used 
in the literature, yet none of them is sufficient to simulate braneworld and we have to develop 
them further. In this thesis, the generalized harmonic formalism will be employed to evolve the 
branewrold spacetime. Thus we put the introduction and the development of the BSSN formalism 
to appendix A. 

1.6 Numerical Methods 

The equations of motion of gravitational theory are non-linearly coupled partial differential 
equations (PDEs). Due to the non-linearity and the complexity, it is not very realistic to study 
the full dynamics in closed form, especially the behaviour at high energy regime where the fields 
are so strong that perturbation methods do not apply. We here use a numerical approach. In this 
section we introduce hnite difference approximation (FDA) methods to solve the PDEs. The focus 
is on the various tests to distinguish numerical solutions from numerical artifacts. 
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1.6.1 Finite Difference Approximation 

To demonstrate the concepts in a less abstract way, let us consider the following model problem, 
which is non-linear wave equation in flat 3 + 1 dimension spacetime under axisymmetry with source 
term / (which does not depend on the wave function <F). This model problem includes a few 
features that are important for numerical calculation in braneworlds. The equation is assumed to 
be (—dtt + d xx + d yy + d zz ) <f> + <f > 2 = / in Cartesian coordinates, or 

+ dp P + —dp + 9 2Z j < F + < L 2 = /, (1-79) 

in cylindrical coordinates (t. p, <f>, z) that are adapted to the axisymmetry. Therefore the axisym¬ 
metry implies d^ = 0, which has been applied in (1.79). Let us assume the spatial domain is 
P G [0, Pmax] i Z G [0, 2max]- 

The whole domain, both spatial and temporal, is divided into discrete grids (or meshes). In 
principle this division can be arbitrary, as long as the grid/mesh elements are small. The meaning 
of “small” is going to be clear by the discussion in section 1.6.2. To be more specific and to honor 
simplicity, here let us employ uniform grid. Therefore the spatial domain can be 

Pi = (i — l)Ap, i = 1,2,... , rip where Ap = /Jmax ; (1.80) 

Tip — 1 

Zj = (j-l)Az, j w 1,2,..., n z where Az = ' m ^ x . (1.81) 

Tl z -L 

For simplicity, let us choose Ap = Az = h. The time domain is also discretized and the time 
interval between two subsequent discretized time levels can be expressed as At. A t/h is called the 
Courant factor. 

We use notation 

= *(t n , Pi ,Zj) = $((n - 1 )A t, (* - l)Ap, (j - 1 )A z), (1.82) 

and similar notation for function /. We replace the differential operators by their FDA operators 
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with second order accuracy: 


dpp* 

-+ 

*?+ij ~ 2*? tj + 

(1.83a) 

h 2 

d P * 

-+ 

A-ij 

(1.83b) 

2 h 

d zz * 

-> 

$ n ■ 11 - 2<F n ■ + , 

1,3+1 t,3 ' i,3~ 1 

(1.83c) 

h 2 

d tt * 

-> 

$"+l _ 0<b n . 4- (J)" -1 

1,3 *,3 1,3 

(1.83d) 

(Ah) 2 


The FDA operators are obtained by Taylor expansions such as 

h 2 


<t> n , . = <f> n . 4- h(f> _<F> „„ -|_$ -|_$ 

'■+1 -J ^*,3 ' n ^,P ' 21 >00 ' gl >PPP w 4j 1 ■/ 


h 3 „ h 4 . 

,ppp ' 


+ 0(h 6 ), 


which yield 


ci>" , . - 2<h n . | <I>" , . h 2 

a +w _ ; -bj - a $ + + o(h 4 ) 

jj2 — U PP ^ + Y2^' PPPP ^ l '' 


(1.84) 


The term jo*, pppp + 0(h 4 ) = 0(h 2 ) is the difference between the exact operator and the FDA 
operator, which is called truncation error. When h is small (so that the truncation error is not 
significant), the differential operators can be replaced by their FDA counter parts. Other FDA 
operators in (1.83) can be obtained similarly. The discretized PDE reads 


. *?+u-2. i 

I 7 0 I 


m 2 


h 2 


2 h 


$i J+ 1 - 2 <b( l ■ + $ 


i,j w 


h 2 


+ ( $ hi ) 2 = fhr ( L85 ) 


Now we are ready to introduce the general notations to make the discussion clearer. A set of 
PDEs, such as equation (1.79), can be collectively denoted as 


Lu = /, 


( 1 . 86 ) 


where L stands for differential operators and all other operations, u stands for the fundamental 
variables (the unknown functions) to solve for, and / stands for the terms in the equations that 
do not include u. In equation (1.79), u = $, and Lu = L$ = [—d tt + d pp + ^ d p + d z ^j <f> + 4> 2 . 
The discrete FDA operators, such as equation (1.84), can be collectively denoted as 

A<1> = £<h + h p ■ E<f>, (1.87) 
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where A stands for the FDA version of the exact operator S. h p means that the approximation 
level is of p-tli order in h. E stands for the error operator — more specifically, h p E<h is the error. 
Using (1.87), we can discretize (1.86) as 


L h u h = f h , (1.88) 

where h is to label resolution. An example of (1.88) is (1.85). 

In (1.85), the approximation is of second order in h. Generally the approximation order of L h 
is p, which can be formally expressed as 


L h = L + h p E. (1.89) 

From the discussion above, one sees that the validity of FDA needs to be built upon the 
following two assumptions: ( 1 ) the funtion $ is smooth; ( 2 ) h is small, so that the truncation error 
is not significant. 

However, these two conditions are not sufficient to guarantee that the numerical result u h is 
actually a approximation of the exact solution u. Therefore systematic test mechanisms need to 
be developed to distinguish numerical solutions from numerical artifacts. 


1.6.2 Tests 


First, often it is neither practical nor necessary to let equation (1.88) be satisfied exactly. 
Instead, (1.88) is considered to be satisfied when residual r h = L h u h — f h is sufficiently small. 
Again, “small” does not have any measurable meaning yet. 

Multiplying equation (1.85) by pi , we get the following equation 


p i ^n+l_ 2 ^n j + ^n- 1 ) Pi ^ +1<J - 2^ .+^_lj) , ^i+U ~ 


(A ft ) 2 


h 2 


2 h 


i• a /',//:,• (i-90) 


(1.85) and (1.90) share exactly the same numerical properties, such as convergence, smoothness, 
regularity, etc. But the two residuals have different numerical values. Therefore, the residual being 
“small”, has no absolute meaning. 

This feature can be expressed in a more abstract way as: Lu = f and g ■ Lu = g ■ f have the 
same numerical properties. Here g is a non-zero, smooth function over the domain. For example g 
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can be an arbitrary non-zero constant to make the residual take any value. Therefore, the absolute 
value of residual does not have any meaning. So how to distinguish between a numerical solution 
and a numerical artifact? And how small is the residual to be considered sufficiently small? These 
questions will be answered by the following analysis. 

Assume the numerical result u h that satisfies L h u h = f h -\- r h is obtained, where r h is the 
residual. Generically, u h is a numerical solution, if the following equation is satisfied when u h is 
substituted back into equation ( 1 . 86 ) 


lim Lu h - f h = 0. (1.91) 

h-> 0 V ' 


Let us see what it means 


Lu h - f h = L h u h - h p Eu h - f h = r h - h p Eu h = r h + 0[h p ). (1.92) 

Therefore (1.91) is satisfied, if r h is negligible compared to h p Eu h (in this sense r h is small). 

However, technically it is impossible to apply a continuous operation L to discrete function 
u h , and then eq. (1.91) can only be understood formally. Instead, u h is considered a numerical 
solution, if 


lim rj 1 = 0 , 
h^O 

where rj 1 = LiU h 


f h , where L 7 ^ L h that satisfies lim = L. 

h— s -0 


(1.93) 


Since L is independent of L h (a different discretization), rj 1 is called independent residual. 

For the model problem, we can use the following discretization as the independent discretized 
operators 

(1.94a) 
(1.94b) 
(1.94c) 
(1.94d) 

This discretization is different from (1.83), and is also of the second order accuracy. 


d r r*h —}■ 


d r $ - 

du$ -»• 


2$ n . - 

r ^i+u 

+ ^?+2J 

— <3) n 

^+3,j 



h 2 


CO 

< 0 . 

- 4^?+i 

,j + ®i+2J 



2 h 


5 

- 

c,d>™ 

0y -i,j + l 

+ 4$ ”j+2 

— <f> n 

P i,j+3 



h 2 


2 - 

- 

f 4$”- 2 - 

CO 

1 

s -*r 


(A hf 
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In general, the approximation order of Lj is denoted as m, therefore 

Li = L + h m E l = L h - h p E + h m Ei, (1.95) 

rf = L^u h -f h = (L h - hPE + h m E 1 )u h - f h = r h - h p Eu h + h m E lU h . (1.96) 

Again, here it is required that 11 r h 11 is negligible compared to min (11 h p Eu h | . h m Eiu h 11), therefore 
the independent residual rj 1 converges to zero at min(p, m)-th order. Here ||rt|| is the norm of u. 

For the model problem, p = m = 2, therefore the independent residual behaves as a second 
order quantity: when h decreases to h/2, the independent residual decreases to r^ 1 ^ 2 ' 1 = \r^. 

Note , the independently discretized operators L{’ can be very different from the discretized 
operators L h used to obtain the solution. L f { and L h do not need to be of the same method. For 
example, one can use finite element method or spectrum method to obtain the solution, but use 
FDA as independent operators to evaluate the independent residual. 

1.6.3 Tests for General Relativity 

For a numerical problem, often there are a certain number of equations to solve, for an equal 
number of fundamental variables (the unknown functions). If the number of equations is less than 
the number of unknown functions, in principle there are no unique solutions. On the other hand, 
in GR, the number of equations is greater than the number of unknown functions. In this case the 
redundant equations are called constraints. 

As an example, in 3 + 1 formalism of GR, there are six functions 7 ij to be solved for, by 
solving the six evolutionary equations. The other four equations are the Hamiltonian constraint 
and momentum constraints. Analytically, if the constraints are satisfied initially, the consistency 
(Bianchi identity) guarantees them to be satisfied at all times, as long as the evolutionary equations 
are satisfied during the evolution. However, numerically there are always small violations to 
the constraints, and there is 110 guarantee the violations are controllable. Therefore, for general 
relativity, the constraints need to be tested as well. i.e. in order to make sure all the components 
of Einstein’s equations are satisfied, both the independent residual test and the convergence test 
for constraints are needed. 

Equivalently , in the case a certain formalism of GR is employed to obtain the numerical results, 
the results can be substituted into another formalism of GR to produce residuals, and the residuals 
should converge at the expected order. For example, one can use generalized harmonic formalism 
to obtain the solution, and then substitute the solution into original Einstein’s equations to get 
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residuals, and check whether the residuals converge as expected. 
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Chapter 2 


Characteristics in the Braneworld 
Spacetime 


The presence of the brane imposes interesting new physics. This chapter is devoted to develop 
the formalisms to describe the following topics associated with the brane. 

Israel’s junction conditions impose cusps in some metric components, which serve as boundary 
conditions for these metric components. In this chapter we will discuss the boundary conditions 
of the remaining metric components. We will also discuss other effects of Israel’s conditions such 
as the smoothness of the apparent horizon across the brane. 

The main goal of our study is to numerically simulate the process of black hole formation. The 
definition of a black hole relies on the global causal structure of the spacetime. We will discuss 
how a black hole can be defined in the braneworld. In the braneworld, the causal structure of 
the spacetime is determined by the spacetime geometry in 5D, therefore the 4D apparent horizon 
and event horizon on the brane should play no direct role in braneworlds. However, since the 4D 
brane is all one can observe, we will study the 4D quantities as well, and compare them with the 
results in GR to see the observable difference from GR. We will also study the relation between 
the horizons on the brane and the horizons in the bulk. 

Energy in GR is not a locally defined quantity since the energy “density” can be of any value 
[94]. However, if the system presents asymptotic translational symmetry in time, in certain cases 
quasi-local energy can be defined, such as ADM energy in asymptotically flat spacetime. A more 
general formalism of energy obtained from Hamilton-Jacobi analysis [95] [87], will be directly 
employed in the braneworld to obtain the total energy of the system. 

There is also energy exchange between the brane and the bulk. In this chapter we will present 
our preliminary study on this topic. 
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2.1 Boundary Conditions at the Brane 

In this section we discuss the properties and the gauge freedom in the boundary conditions at 
the brane. 

The vacuum solution of the braneworld is eq. (1.38) 

ds 2 = e - 2M/i (V dx^dx") + Ay 2 , (2.1) 

where y G (— 00 , 00 ), and the brane is located at y = 0. For the general case (non-vacuum), we 
setup the coordinate system (x a , y), where y is the extra dimension. Latin indices (a, b ,...) are for 
the coordinates on y = constant surfaces, and their values take 0,1, 2,3. Greek indices (/.i, v ,...) 
take 0,1,2,3,4, and are used for the coordinates of the whole spacetime. Therefore the metric is 

ds 2 = g a bAx a dx b + 2g ay dx a dy + g vy dydy. 

The coordinate y is set to inherit the following properties: y = 0 is where the brane is located, 


and y is adapted according to the Z 2 symmetry, such that the metric components: 

9 ab(x a , -y) = g a b{x a , y), ( 2 . 2 ) 

9yy( x >—y) = 9yy( x (2-3) 

and g ay {x a ,-y) = -g ay (x a ,y). (2.4) 

Under this coordinate choice, Israel’s first junction condition is simply g a b\ y= o+ = 9ab | y=0 - = hab. 


where h a b is the intrinsic metric on the brane (expressed in the coordinates x a ), induced from 
the bulk metric on either sides of the brane. Israel’s second junction condition imposes conditions 
on the extrinsic curvature K. a b (of the brane embedding in the bulk). These conditions can be 
translated into conditions on dg a b/dy, which will serve as the boundary conditions for g a b- 

g yy and g ay are not related to Israel’s conditions. Rather, since the braneworld spacetimes are 
“one-sided” (see Sec. 1.3.2 and Sec. 1.3.4), in general there is 110 generic reasons to require the 
y-coordinate lines (the intersection of the x a = constant surfaces) to be perpendicular with the 
y = 0 surface (the brane), which means g ay \ y _ 0+ 7 ^ 0. Taking y —> 0 in (2.3) and (2.4), we get 



to 

II 

O 

1 

II 

to 

II 

0 

+ 

(2.5) 

and 

9ay y=0 - ~ 9ay y=Q+ ■ 

(2.6) 


37 




2.2. Apparent Horizon 


However, since g ay | Q+ 0 in general, it means g ay | Q are not defined. This is because only the 
induced metric on the brane and the extrinsic curvature of the brane are important, while g ay are 
not needed in defining the induced metric and the extrinsic curvature of the brane. Similarly, g yy 
is not needed in defining the brane geometry either, which means generically g yy is not defined, 
although it could be defined as g yy \ Q = g yy | Q _ = g yy \ 0+ - We have 


generically, g yy are not defined on the brane. 

Although generally g^ y are not defined, it is convenient to choose the coordinates at the brane such 
that the y-coordinate lines are perpendicular with the brane. We call this gauge as perpendicular 
gauge. Under this coordinate gauge, we can then define 

9ay\y =0 = 0. (2.7) 

This coordinate gauge has desirable properties such as the smoothness of apparent horizon that is 
going to be in Sec. 2.2. 

Note that g yy \ 0 is still unconstrained. 


2.2 Apparent Horizon 

An apparent horizon is needed to monitor the evolution of spacetime if a black hole is present 
during the evolution. This section is devoted to studies on apparent horizons in braneworlds. 

2.2.1 The Definition 

The definition of apparent horizon can be found in standard texts [34, 35, 37]. Let n a be future 
directed timelike unit vector normal to t = constant hypersurfaces. Let S denote a closed (d — 2) 
dimensional spatial surface within a t = constant hypersurface, and s a is its unit normal vector 
pointing towards the outgoing direction, which is within the same t = constant hypersurface. The 
induced metric on S is then (not to be confused with the m a defined in (1.3)) 

m aP - g ap +n <y n P _ s a s f}_ ( 2 . 8 ) 
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Let v a be an arbitrary vector field, the relative change rate in the area elements of S along v a is 


[37] 




— \/rri — mT /L v n~ini/ — mT VTrv, 
\Jm 2 


(2.9) 


where m is the determinant of . The final expression is the expansion of v a along S. 

For the (d — 2) dimensional spacelike surface S, there exist two null curves orthogonal to this 
surface. Let us denote the two future directed null vectors tangent to these two null curves as ^la- 
The relative change rate of the area elements of S along the null geodesics congruences produced 
by % are then 

e± = m a,3 v Q %. (2.io) 


A trapped surface is an S whose 0± < 0, which means the null geodesics congruences produced 
by both + ln and ~l M drag S to the surfaces with smaller area elements. An S with 


9 + = 0 


( 2 . 11 ) 


is called a marginally outer trapped horizon (MOTH), whose area elements stay the same under 
the Lie-dragging of + l M . The MOTH is not unique in a given spacetime since there can be other 
MOTHs within a MOTH. An apparent horizon is defined as the outermost MOTH. 

Now let us construct ^la- Because timelike normal vector n a is orthogonal to the spacelike 
normal vector s a , the two future directed null vectors orthogonal to S are 

± l a = n a ±s a , ( 2 . 12 ) 

where + l a is outgoing and ~l a is ingoing. Substituting this into (2.10), we obtain 

0± = ( np ± sp) = ±D a s a — I\ + s a s l3 K a p, (2-13) 


where D a denotes the covariant derivative associated with 7 a p, and K a p is the extrinsic curvature 
of the t = constant hypersurface. Therefore, the follow equation is satisfied at the apparent horizon 

0+ = D a s a — K + s a s /3 K a p = 0. (2-14) 

The above discussion shows that the definition of apparent horizon relies on the choice of 
t = constant hypersurfaces. For a given spacetime, different slicing conditions can result in drasti¬ 
cally different apparent horizons. Taking Schwarzschild spacetime as an example. In Schwarzschild 
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coordinates or the isotropic coordinates, the apparent horizon coincides with the event hori¬ 
zon, where the physical singularity is inside of the apparent horizon. On the other hand, the 
Schwarzschild spacetime can be sliced in such a way that there is no apparent horizon [33]. 

2.2.2 Smoothness of the Horizon 

Israel’s conditions impose cusps to some components of the metric. In this section we will study 
whether this affects the smoothness of apparent horizon. The smoothness of an apparent horizon 
can be studied via s“ by asking whether s° is continuous across the brane. Rewriting 0 + = 0 as 


D a s a = I< - s a s f} K a/3 , (2.15) 

we then apply Gauss’ theorem (in curved space) on the t = constant hypersurface. By the same 
procedure to derive the junction condition of electric field across a surface 6 , we can find the junc¬ 
tion condition for s“ across the brane. If the right hand side of (2.15) is finite (by the Z -2 symmetry 
with respect to the brane, this condition is true if the t = constant hypersurfaces are chosen to 
be perpendicular with the brane), then the integration over an infinitesimal layer across the brane 
vanishes. Therefore the component of s a that is perpendicular to the brane, is continous across 
the brane. This continuity of the perpendicular component, together with the Z -2 symmetry with 
respect to the brane, imply that the perpendicular component of s a must be zero. Therefore s“ is 
continuous, i.e. 


the direction of an apparent horizon is continuous across the brane, if the slicing condition is 
such that the t = constant hypersurfaces are perpendicular to the brane. 


Under the coordinates setting in Sec. 2.1, the slicing condition is expressed as gt y = 0. In 
particular, the perpendicular gauge (2.7) satisfies the slicing condition. 

2.2.3 Apparent Horizon on the Brane and in the Bulk 

Generically the causal structure is determined by 5D geometry. However, only the brane 
quantities are directly observable, therefore we study the relation between 4D and 5D quantities. 

6 The procedure to derive E_ — E_ oc a from V • E oc (>. where p is volume charge density and a is areal charge 
density of the singular layer. E is the electric field, and E+ is the electric field on one side of the singular layer, and 
E _ is the electric field on the other side. 
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The question we try to address in this subsection is whether the apparent horizon seen on the brane 
(which is calculated based on the brane geometry only), and the bulk apparent horizon (which is 
calculated based on the bulk geometry), agree with each other on the brane. This can be studied 
via the expansions of the outgoing null geodesics congruences on the brane and in the bulk 

©brane = + < r W* < r W* - ( V) V a (<% + %) 

= < r W^V Q + (r) S/3 ) (a) 

= (r) m a/3 Vo (np + sp) (b) 

= (m a/3 - n Q n / 3 )V Q (np + sp) 

= ©buik-n“n 0 y a {np + sp). (2-16) 

where m= h a & + * r Vi“ is the projection operator that projects to the (cl — 3)- 

surface on the brane, and D is the covariant derivative associated with the brane metric h a p. 
Anything with a superscript (r) is a quantity defined only on the brane. The vector n“ is the 
unit normal vector that is perpendicular to the brane. Assuming perpendicular gauge (2.7), we 
have sp = C ) S p and np = ( r )np on the brane, which are used in deriving eq. (b) from eq. (a). 
The difference between the two 0’s, even at the apparent horizon where ©bulk = 0, is generally 
non-zero. i.e. generally these two apparent horizons do not agree. Therefore, we will study the 
relation between event horizons in the next section. 


2.3 Event Horizon 

2.3.1 Event Horizon in the Braneworld 

In this subsection we examine whether event horizon is well-defined in the spacetime of the 
RSII braneworld, and discuss how to define an evolution problem. 

The global causal structure of AdS spacetime and its Poincare patch was introduced in Sec. 1.4. 
The spacetime background of RSII braneworld is to take the z > 1 portion of the Poincare patch. 
Therefore the z = 0 boundary is eliminated from the Penrose diagram, and the global causal 
structure is the same as that of Minkowski spacetime. In particular, similar to Minkowski space¬ 
time, there is no signal travelling to spatial infinity and then coming back within a finite local 
proper time, and Cauchy surfaces exist [6]. To discuss the Cauchy surfaces, we define the future/- 
past horizons as (see Sec. 1.4.4). Considering a t = constant hypersurface 
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(which is spacelike), any future causal curves coming from the past horizon will either go across 
this hypersurface, or hit the brane and then are reflected to travel back into the bulk (due to 
the Z 2 symmetry with respect to the brane) which eventually go across this hypersurface. Simi¬ 
larly, all the past causal curves of the future horizon go across this hypersurface. Therefore, any 
t = constant hypersurface is a Cauchy surface, because all the developments of any past event are 
captured by the surface, and all the future events can be predicted by the data on this surface. Or 
more precisely, all the inextendible future causal curves of the past horizon and all the inextendible 
past causal curves of the future horizon, go across the Cauchy surface [32]. Therefore, a Cauchy 
problem (an initial value problem) is well-defined. 

The z > 1 branch has the horizons which are not true infinities. Accordingly, the definition 
of event horizon is modified as below. In fact, the definition of event horizon in asymptotically 
Minkowski spacetime can be directly carried over, while the only change is to replace the notion 
of null infinities in the definitions by the future/past Poincare horizons defined by eq. (1.67). The 
event horizon is the boundary of the spacetime region which can not be connected to the external 
world by future oriented null geodesics, i.e. the event horizon is the collection of “the points of 
no return”. The “external world” needs to be defined. Similar to the case of asympototically 
Minknowski spacetime, the external world is defined as the past of the future Poincare horizon, 
therefore the event horizon is the boundary of the spacetime separating the region that can be 
connected to the future Poincare horizon by future null geodesics from the region that can not. 
The future Poincare horizon, on the other hand, are defined as the “infinities” (as measured by 
the Killing parameter t) of future null curves departing from the external world, i.e. there is a 
circular argument in these definitions. The circle can be ended by physically identifying certain 
spacetime region(s) as the external world. Therefore, as long as the external world can be physically 
identified, the event horizon can be defined, and this argument applies to any spacetime (i.e. it is 
not limited to the case of the Poincare patch of the AdS spacetime). 

2.3.2 Event Horizon on the Brane 

The event horizon is defined in Sec. 2.3.1. In the following we will study whether the event 
horizon based on the brane geometry is well-defined. The causal structure of the braneworld is 
determined by the 5D geometry, rather than the 4D geometry of the brane. The null geodesics gen¬ 
erated by the outgoing, future oriented null vectors on the bulk event horizon, form the boundary 
(the event horizon) separating the spacetime region that can be connected to the future Poincare 
horizon by future oriented null geodesics, from the region that can not. Let us consider the out- 
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going future oriented null vectors within the brane originated from the intersect of the bulk event 
horizon and the brane. If the future oriented null geodesics generated by these null vectors will 
stay on the brane forever, then these null geodesics are the boundaries separating the spacetime 
regions that can be connected (by future oriented null geodesics) with the future Poincare hori¬ 
zon from the regions that can not. i.e. (1) these null geodesics stay on the brane forever; (2) 
these null geodesics are the boundaries separating the spacetime regions that can arrive the future 
Poincare horizon from the regions that can not. i.e. these null geodesics form the event horizon 
on the brane. Therefore, the key for the well-definedness of the event horizon on the brane, is 
to study the relation between the null geodesics produced based on the brane geometry and the 
null geodesics produced based on the bulk geometry, generated from the same null vector which 
initially lies within the brane. This relation can be described by the extrinsic curvature of the 
brane. 

2.3.3 Extrinsic Curvature as Geodesics Deviation 

This subsection (Sec. 2.3.3) is the foundation of the study on the relation of event horizons. 
The work in this subsection is a “re-discovery” of the Gaussian curvature (see, e.g. [80]) and the 
Gauss-Weingarten equation [81] in differential geometry. 

The motivation is as follows. Let us examine how to measure the extrinsic nature of the 
embedding of a hypersurface E into a higher dimensional space M. E is considered flatly embedded 
into M, if E appears to be flat in M, in the sense that an arbitrary straight line as seen in E is 
also a straight line as seen in M. Since straight lines are geodesics, the above statement can be 
more precisely rephrased as: if the hypersurface is flatly embedded, an arbitrary geodesics in E 
(consistent with the hypersurface metric "/ap) will also be a geodesics in M (consistent with the 
bulk metric g a p)- For non-flat embedding, these two types of geodesics are not the same in general 
'. Therefore, this motivates us to describe the extrinsic curvature as the deviation of the geodesics 
defined in E from the geodesics defined in M for a vector initially lying on E. This point of view 
to describe the extrinsic curvature, is referred as geodesics point of view (GEP for short). 

On the other hand, the embedding described via the covariant derivative of the unit normal 
vector along E, as what has been done in eq. (1.5), is referred as normal vector point of view (NVP 
for short). 

Generally if a d — C dimensional sub-manifold E is embedded in a d dimensional manifold M, 

7 As an example, we can think of a sphere S 2 : x 2 + [y — l) 2 + z 2 = 1 embedded in R 3 . Take a vector lying on 
S 2 : (d x y at (0,0,0). The geodesics produced by it on S 2 is the equator, while the geodesics produced by it in R 3 
is the x-a.xis. 
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C is called the co-dimension. In this subsection, we will prove NVP and GEP are equivalent in 
the (7 = 1 case, in the sense that the embedding studied from GEP (shown below) will also lead 
to the same definition of the extrinsic curvature as eq. (1.5), which is defined from NVP. In the 
RSII braneworld, there is only one extra dimension, thus the co-dimension is C = 1. For C > 1 
case, please refer to appendix B. 

The basic idea is to study the two geodesics generated by an arbitrary vector T a G E, via 
equations T l3 DpT a = 0 and T^V/gT" = 0, where D is the covariant derivative in E and V is the 
covariant derivative in M. Then we compare the two geodesics to get the difference, which can 
describe the embedding nature of E in M. We adopt the following approach: rewrite T 1 ^ DpT a = 0 
as T^V f}T a = leftover, then the leftover is the deviation of the two geodesics. Let us denote the 
unit normal vector of E as n a whose length square is e = n a n a = ±1 where e = 1 if the extra 
dimension is spacelike and e = — 1 if the extra dimension is timelike. The reduced metric of the 
hypersurface is 

7 a/3 = g al 3 - e n a np. (2-17) 


For a general tensor T 01 '" 0 ^ ^ G E, the covariant derivative associated with the metric 7 a ^ 
is [34, 36] 

1 .../ 9 i = 7 \- 7 / 3 i £ ‘ 7 /V 1 ,T 5 -\... £i . (2.18) 


7T) rpOL\...Otk 

fi\. 


Therefore V T“ 6 E, we have 


D a T» = 7/7 wr, 


(2.19) 


which tells us that the geodesics generated by T a in E is just 


0 = T a D a T» = 


7 M 7 T“V a T 7 , 


( 2 . 20 ) 


or 

7 ^ v T a S7 a T v = 0. (2.21) 

On the other hand, we have the following identity which can be obtained by the fact that n a T a = 0 

n^n u T a V aT 1 ' = -n^T a T v V a n v . ( 2 . 22 ) 

From (2.21) and (2.22), we obtain 

T“V„T M = ( 7 ^ + e n»n v ) T Q V a T v = -en^T a T^V a n p . (2.23) 
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We can now define the deviation of the two types of geodesics equation as the right hand side of 
equation (2.23). It is then clear that the deviation vanishes, if and only if 


T a T^V a np = 0. (2.24) 

This result is for arbitrary T a defined on E, and only the contraction with T a appears in this 
expression, which means we can use the following quantity to describe the embedding of E in M 

K al 3 = -7/7/V^n,. (2.25) 

i.e. NVP and GEP are using the same quantity ( K a p ) to describe the embedding. The deviation 
equation (2.23) in terms of the extrinsic curvature is now rewritten as 


rV a F = e n^T a T p K aP , 


(2.26) 


from which one see that the deviation is in the perpendicular direction (n M direction). 

2.3.4 The Relation between the Event Horizons 

The work in this subsection (Sec. 2.3.4) was first independently carried out in [81]. 

For the braneworld, the 3 + 1 brane is embedded into the 4 + 1 dimensional bulk. Using the 
notations in braneworld, the geodesics deviation equation (2.26) is now 

/V a / = jT+ a + /3 /C a/ 3, (2.27) 

where v a is the tangent vector of an arbitrary geodesics within the brane. n is the unit normal 
vector of the brane, and e = = 1 has been applied. 

To study the event horizon relations, we focus on the case where v a is a null vector. In RSII, 
K. a p is related to brane content by Israel’s junction condition (1.29) 

= 2 kd ^ + T a p - h a p 2 ^ > (2.28) 

which implies the deviation of null geodesics is 

v a V a v^ = ^k d n^v a v P T a p. (2.29) 
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Therefore the deviation of the two geodesics amounts to whether v a v l3 T a /3 vanishes. Generally 
the right hand side of eq. (2.29) does not vanish (otherwise it is a new energy condition for the 
matter. The discussion of energy condition is beyond the scope of our project). However, in case 
Tfxv = 0, when the matter on the brane vanishes, the two geodesics coincide. In another word, 
when the matter is absent, at the intersect of the bulk event horizon and the brane, the future 
oriented null geodesics that are produced by the null vectors lying in E, will stay on the brane 
forever. Since these null geodesics are on the event horizon of the bulk, they are the boundaries 
separating the spacetime regions that can be connected (by future oriented null geodesics) with 
the future Poincare horizon from the regions that can not. i.e. (1) these null geodesics stay on the 
brane forever; (2) these null geodesics are the boundaries separating the spacetime regions that 
can arrive the future Poincare horizon from the regions that can not. Therefore, from the brane 
point of view, they form the event horizon on the brane. i.e. 

an event horizon on the brane is well-defined when there is no matter around the horizon. 

For gravitational collapse processes, if the systems eventually reach the stationary states that 
the matters either go into the black holes, or get bounced back and travel towards spatial infinity, 
then there are no matters at the horizons and the event horizons on the brane are well-defined. 

2.4 Energy in the Braneworld 

In order to quantitatively describe the spacetime evolution, and gravitational interaction be¬ 
tween the brane and the bulk, we need to introduce certain quantities, such as energy. However, 
there is no local definition of energy in general relativity. Instead, there have been many attempts 
to define quasilocal energy in general relativity, and many of these definitions are only well-defined 
in a certain background. In this section, the definition developed by Brown and York obtained 
from a Hamiltonian-Jacobi analysis of the gravitational action [95] (Hawking and Horowitz also 
gave a similar derivation [87]), is applied to the braneworld. 

2.4.1 Total Energy 

In this subsection, we introduce the energy defined in [87, 95]. In the next subsection, we will 
apply this definition to the braneworld spacetime. 


46 







2.4. Energy in the Braneworld 


In general the action of a gravitational system is 


1(9,*) 


l \— 

Jm L 167t Gd 


+ -^m(g, *) 


1 

87t Gd 



(2.30) 


where M is the spacetime manifold, Jff is the extrinsic curvature of the boundary dM embedding 
into M. «5? m is the Lagrangian for all the matter fields and the matter fields are collectively denoted 
as f I>. 

Let us choose the boundary DM = where BQ = Ute[ii,t 2 ] S t as shown in 

Fig. 2.1. Sf is the closed (d — 2) dimensional spatial surface embedded in each E t , and Q is the 
single parameter to characterize the family of the enclosed (d— 2) dimensional surfaces in each Et. 
When Q —> oo, Sf goes to the spatial infinity boundaries. B® is setup such that its normal unit 
vector is perpendicular to n M (so that n M lies within B®). 



Figure 2.1: Manifold M and its boundaries. The boundary dM is composed of E tl , E ta , and B®. 
On the diagram, the boundaries are shown to be the upper/lower surface of the cylinder, and the 
side surface of the cylinder. Here B® = |J tg r ti t 2 ] > where S® is the enclosed (d — 2) dimensional 

spatial surface embedded in each E t , and Q is the single parameter to characterize the family of 
the closed surfaces in each E 4 . Choose B® such that its normal unit vector Q^ is perpendicular to 
(so that lies in B®). 


The Hamiltonian is then [87, 95] 


H = 


[ (an + $ V M V ) + -±— QA K ^ V ~ + « • {d ~M , (2.31) 

7s t oirGd Js , Q L 


where e = n M n p = — 1. ( d 2 \ a p is the extrinsic curvature of Sf embedded in E t , and ( d 2 \ is its 
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trace. H and M v are Hamiltonian constraint and momentum constraint defined as 


U = TGttCTi (~ K2 + + e • {d ~ 1)R + 16 * G dp) , (2.32) 

M v = -^—DJKy^ - K^) - S^. (2.33) 

onGd 


The Hamiltonian constraint H and momentum constraints M are implied by Einstein’s equations, 
therefore vanish for physical configurations, and should be dropped in (2.31). 

Hamiltonian (2.31) diverges in general. However, it is the physical Hamiltonian that matters 
[87]. The physical Hamiltonian is H — H, where H is the Hamiltonian of a background spacetime. 
We denote background quantities by a bar (“). The background is a static solution, then its 
contribution is only 


1 

87t G d 



G ~ 2 \ 


where the integration is over a closed surface in the background spacetime that is isometric to the 
Sf (that is the closed surface chosen in the physical spacetime). 

The physical energy is the physical Hamiltonian 


E = H — H = lim <L \-e QJKy ^ - K^ v )p v + a ( (d ~ 2 \ - G ~ 2 \\ 1 , (2.34) 

Q— >oo 87rGd Js® V / J 

where the Hamiltonian constraint H and momentum constraint A4 are dropped, since they vanish 
at physical configurations. Also, at spatial infinities where there is asymptotic time translational 
symmetry, the lapse function and the shift functions go to the form of the background spacetime. 

However, as proved by Shi-Tam [90], the definition does not yield a positive definite energy 
except for the time symmetric case (A' ; ,„ = 0). Also, the definition is gauge dependent. Yau and 
Liu [91-93] defined another formula for energy, which is gauge independent, and can be positive 
definite under certain conditions. However, our initial data is time symmetric which means the 
definition by Brown-York and Hawking-Horowitz is sufficient. Also, we can still use the definition 
by Brown-York and Hawking-Horowitz during evolution since the energy is characterized by the 
asymptotic behaviour at spatial infinities, which is not affected by finite time evolution (i.e., the 
local behaviour is not able to propagate to spatial infinities within a finite time evolution). 
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2.4.2 Total Energy in the Braneworld with Axisynnnetry 

We assume the energy (2.34) is well-defined in the braneworld. In the braneworld, the spacetime 
background is 

ffi- / _ _ _ \ 

ds 2 = — ( — dt 2 + dr 2 + dz 2 + r 2 (d# 2 + sin 2 9d(j > 2 ) J, where z> £. (2.35) 

A bar is used to denote the quantities associated with the background. For this background, the 
lapse a = l/z and the shift /3j, = 0. Therefore, (2.34) reduces to 

E =« li ?»sk£ P a ( w " >t -""*)■ < 236 > 

which will be the definition of the energy of the braneworld. 

To calculate the energy, we need to set a closed surface family S® that goes to spatial infinity 
as Q —> oo. The two requirements on defining the family are: (i) S® goes to spatial infinity as 
Q —* oo, and (ii) the closed surface is smooth to a certain degree so that the extrinsic curvature 
(d ~ 2) k afi is well-defined at any point on S®. In our case where the system has axisymmetry 
(spherical symmetry on the brane) with coordinates (t, r, 9 , <p, z) adapted to the symmetry, we may 
choose the closed surface as, for example, Fig. 2.2(a). Quantities Q , u and v are the parameters 
for defining the closed surfaces. Please refer to the capture of the figure for the details. 

Without loss of generality, the spatial metric of any t = constant slice can be 

d l 2 = ^ [e 2 - 4+2S (dr 2 + d^ 2 ) + e 2A ~ 2B r 2 (dfl 2 + sin 2 9d(j> 2 )] . (2.37) 

The brane is located at z = l. This is the most general spatial metric for the axisymmetric 
configuration because, by taking the symmetry into account, the most general form can take 

Z 2 

d l' 2 = — [?) rr dr 2 + rj zz dz 2 + 2fj rz drdz + r 2 fjgg (d# 2 + sin 2 #d<(> 2 )] . (2.38) 

For a given t = constant slice, everything only depends on r and z, therefore fj rr dr 2 + fj zz dz 2 + 
2fi rz drdz can transform into a conformally flat form. Lastly, the freedom in the conformal function, 
can be used to fix the brane at z = £ [14, 47]. 

To calculate the energy, we embed Sf into the background spacetime (2.35). In general, it is 
not guaranteed that such embedding is possible for an arbitrarily chosen closed surface, although it 
turns out all the closed surfaces considerred in this thesis could be embedded into the background 
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(a) (b) 

Figure 2.2: The closed surfaces to calculate total energy and energy in the bulk. Diagram (a) is 
used to calculate the total energy, z = 1 is where the brane is located. The system is spherically 
symmetric on the brane and r is the radius of the spherical coordinate. The closed surface is 
composed of three segments: segment (I) starts at r = Q and its coordinate length (measured by 
the coordinate z) is v ■ Q. The length of segment (III) is u ■ Q. Here u G (0,1) and v G (0,1) and 
their values are fixed for a specific closed surface family. Segment (II) is an arc (a quarter of a circle 
whose radius is (1 — u) ■ Q)) to connect these two segments smoothly. The small segment below the 
brane is to show that there is another part below the brane which is related to the said part by 
symmetry. Note, the closed surface should goes smoothly across the brane (i.e. the closed surface 
is perpendicular to the brane). In general case this perpendicular surface is not represented by 
r = constant. When the coordinate gauge at the brane is taken to be perpendicular gauge (2.7) 
so that g rz \ _. = 0, this surface is represented by r = constant. Diagram (b) is to add two dashed 
lines (along the brane) onto diagram (a). The use of (b) is going to be explained in Sec. 2.6. 
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spacetime. The embedding is a mapping from Sf in the physical spacetime to its image in the 
background spacetime by keeping the intrinsic geometry of the closed surface, i.e. the intrinsic 
geometry of the image of Sf in the background spacetime is the same as the intrinsic geometry 
of Sf in the physical spacetime. There is a freedom in this mapping (which is “where we put the 
image”), and we fix this freedom by naturally mapping the intersection of the closed surface with 
the brane in the physical spacetime to the brane of the background spacetime 


Z\z=t = t 


(2.39) 


For a closed surface in the physical space with metric (2.38), the embedding into the background 




(a) (b) 


Figure 2.3: The embedding of a closed surface. Fig. (a) is the closed surface in the physical 
spacetime (2.38), and Fig. (b) shows its embedding into the background spacetime (2.35). The 
freedom of the embedding is fixed by mapping A to A. 


space (2.35) is demonstrated by Fig. 2.3. The embedding boils down into the following two condi¬ 
tions: 

(1) On the r — z plane, an point B (see Fig. 2.3) on the closed surface, represents a 2-sphere 
expanded by coordinates ( 6 , <j >), whose proper area is 47r r 2 fjee/z 2 . The area of the 2-sphere at 
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B (the image of B in the background spacetime) is 4irr 2 / z 2 , therefore 


—gee = —■ 

Z z Z z 


(2) C is a point on the closed surface that is infinitesimally close to B. Let us denote the coor¬ 
dinates of B by (r,z), and denote the coordinates of C by (r + dr, z + dz). The square of 
the length of the infinitesimal line is ( fj rr dr 2 + fj zz dz 2 + 2fj rz drdz) /z 2 . Correspondingly, the 
coordinates of their images ( B and C) are (f,z) and (f + d f,z + d z), and the square of the 
length of the infinitesimal line is (df 2 + dI 2 ) jz 2 . The equality of these two lengths reads 

■yj (f)rrdr 2 + f/ zz dz 2 + 2fj rz drdz) = (df 2 + dz 2 ) . (2-41) 

Here we emphasize that the infinitesimal line is taken along the closed surfaces. The condition 
“along the closed surface” defines how dr and dz are related (also defines how df and dz are 
related). 

The freedom related to “where to put the images” is fixed by mapping A (the point of the closed 
surface on the brane) to A (the point of the closed surface on the brane). 

For metric (2.37), these two conditions reduce to 

2 -2 
ryZj 

—e 2j4-2B = ^ (2.42) 


_L e 2A + 2B ( dr 2 + dz 2 J = J_ ( df 2 + d ^2) 


For background metric (2.35), the lapse funtion to evaluate (2.34) is a = ifz and the shift 
function is /3 M = 0. <f> S Q (a ■ ! d-2 )/c) diverges and the divergence is cancelled by the divergence in 
§ s q (a • ' 2 *fc). Noticing the l 2 /z 2 factor, we examine whether the contribution from segment (II) 

and segment (III) in Fig. 2.2(a), vanishes as Q —> oc. In fact, one can show by direct calculation 
that <j> s Q (a ■ ( d-2 )fc) from segment (II) (the arc) converges to zero as 1 /Q, as long as A and B con¬ 
verge to zero (at any rate) as Q —> oo, so does <f s Q (a • ( d-2 )fc). Therefore <j> s Q a ((d-2) fc _ (d-2)£) 5 
the contribution to the total energy, converge to zero at least as fast as 1/Q. The same argument 
also applies to segment (III). Therefore, only the r = Q segment contributes to the total energy. Or 


within the bulk, only the region r>z, contributes to the total energy. 
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Now we complete the proof by proving the claim “§ s q (a • ( d ~ 2 \) from segment (II) and seg¬ 
ment (III) converges to zero as l/Q”. The calculation on segment (III) is similar to the calculation 
on segment (II). Furthermore, acually segment (III) does not need to exist to define the closed 
surface, if we take u = 0. For conciseness, here we only present the calculation on segment (II) 
(the arc). 

The metric is (2.37), where r and 0 along the arc can be parametrized as r = uQ + pcosy and 
z = 1 + vQ + psiny, where p = (1 — u)Q, x £ [0,7r/2]. A direct calculation gives 

lim I (a ■ = 87t lim / dy e 2A ~ 2B \ (£i 4>&) 

Q—>oo JgQ V / Q—>oo Jq Z^ 

/•"■/ 2 1 

= 8tt lim / dy — (Ci + 6) ■ (2-44) 

< 2 ^°° Jo z 

Here = — r 2 — 2rp cos y + (3r 2 psiny)/z and it is easy to numerically check 



oc 1 IQ. 


£2 = r 2 p (— 3A p + B p ), and A iP and B p converge to zero at least as fast as 1/p since this cor¬ 
responds to A ~ logp and B ~ logp, whilst the reality is A —> 0 and B —> 0 when p —> oc. 
Therefore 


r/2 


dy (6/^ 3 ) 


JO 


converges to zero at least as fast as 


rsj 



dy (r 2 /z 3 ) oc 1/Q. 


This completes the proof. 

Another potential contribution to the total energy is the “cusp” of the closed surface across the 
brane, if there is any. Within the physical spacetime, we can properly choose the closed surface 
such that it crosses the brane smoothly. Under the coordinate gauge g rz = 0, the smoothness 

condition is equivalent to d?’/dz|_ ( = 0, where the derivative is taken along the trajectory of the 

closed surface in the r — z plane. This condition is satisfied by the surface specified by r = Q. 
Within the background spacetime, the smoothness condition is equivalent to dr/d 2 |__ t , = 0. If we 
use z as the parameter of the trajectory in the f — z plane of the background spacetime, we have 


dr/d 2 


dr/dz 

dz/Az 1 
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where the derivatives are taken along the trajectory of the closed surface in the f — z plane of 
the background spacetime. Therefore the condition dr/df|A ^ = 0 is just df/d;z| = _^ = 0. If this 
condition is met, there is no cusp (therefore no contribution to the total energy) across the brane. 
In fact, since we will take Q — * oc eventually, it is sufficient to examine whether this condition 
is satisfied by only keeping the lowest order of A and B in the expressions. However, condition 
df/dz\-_ ( = 0 is generally not gauranteed. A direct calculation from eq. (2.42) and eq. (2.43) 
yields 

df/dz| « r ■ {A + A * + B - B,~) \ z=( , (2.45) 

where « means this relation is true up to the lowest order in A and B. In general, it is hard to 
say whether this term vanishes. However, in the conformally flat space that we are going to study 
in the next subsection (Sec. 2.4.3), the smoothness condition is indeed satisfied (see eq. (2.53)). 
Furthermore, we will obtain a relation between the total energy of the whole braneworld, and the 
ADM energy calculated based on the brane geometry. 

2.4.3 Total Energy in Conformally Flat Space of the Braneworld 

A general way to realize the embedding is to encode the two conditions (eq. (2.42) and (2.43)) 
into numerical methods. However, if we take into account the asymptotic behaviours of A and B 
in the following special case, the embedding and the calculation of total energy can be obtained 
analytically without using numerical methods. 

We consider a special case by assuming B -C A in r z region 8 . In this region, the spatial 
metric reduces to 

dl 2 = ( —e 2A (dr 2 + d^ 2 + dd 2 + sin 2 ddcj) 2 ) , (2.46) 

for which the Hamiltonian constraint is (under the unit t = 87rG5 = 1, which implies 8 ttG 4 = 1 
via (1.37)) 

(d rr + d zz ) A + 4(1 — e 2A ) + ^ ^ + (A tZ ) 2 + (A r ) 2 = 0, (2.47) 

and Israel’s junction condition is 

A, z \ z=l = l~e A . (2.48) 

In the r z region (which is the only region that contributes to the total energy as proved above), 

8 This condition was observed (i.e. a posterior result, rather than a prior assumption) in static star solutions [47]. 
This condition holds in a solution for small static BHs [14]. This condition is also compatible with Hamiltonian 
constraint [48]. This condition will also turn out to hold in one version of our initial data (see Chap. 4). 
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eq. (2.47) reduces to the following by taking only the linear term in A 


(d rr + d zz ) A — 


4 A + 2A : , 



(2.49) 


with linearized Israel’s condition as 



(2.50) 


The Hamiltonian equation (2.49) subject to boundary condition (2.50) can then have an solution 
with closed form 

, (2.51) 

rz 

here we use ~ to emphasize that this solution is only true at r > 2 region. 

Remark: this solution satisfies the linearized Hamiltonian constraint subject to the linearized 
Israel’s boundary condition, and the boundary condition at r —> oo. However, there is no proof 
regarding uniqueness. Therefore, we should justify (2.51) is indeed a solution via the result obtained 
by numerical calculation, which is going to be carried out in Sec. 4.3.3. 

In physical spacetime, the coordinates of the r = Q surface are now (z, #,</>). We will find 
the metric of the corresponding surface in the background spacetime expressed under coordinates 
(z, 9 , <j>) as well. The two conditions of the embedding are now 


r 

z 


( f ') 2 + ( z ') 2 

~2 



(2.52a) 

(2.52b) 


where r = Q is a fixed value, and 1 denotes d/dz. Eq. (2.52) is the ODE set that defines f{z) 

and z(z) for fixed r, subject to the “initial” value z |_ x = 1. By construction, the r = Q surface 

in physical space and r = Q surface in background space have the same intrinsic metric, and the 
coordinates of the surfaces are the same one: ( z,9,(f >). It is then straightforward to carry out 
(2.34). Note the shift functions are asymptotically zeros and the lapse function is asymptotically 
1 /z for background (2.35). 

The solution of the ODE set by keeping the lowest order in A is 


z = z ■exp 


on (z - 1) 
rz 


r = r■ exp 



(2.53) 


It is then a routine work to carry out the calculation of total energy using (2.34). Here we only 
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list some important intermediate results 


Aftotal 




(2.54) 


The factor 2 in the second line is due to the Z-j symmetry with respect to the brane in RSII; 
47T is from the integration over Q is the determinant of the intrinsic metric of the surface: 

Q = e 3A r 2 /z 3 ; and we have used ^ d ~ 2 \ = —2 z/r + 5ai/r 2 , = —2 z/r + 2a\ /r 2 ; r = Q: 

a = e A I z. 

Also, from the expressions one see that, the result would have diverged as ~ Q, if the background 
term ^ d ~ 2 ^k is absent. 


2.5 ADM “Mass” and Hawking “Mass” of the Brane 

Due to the equivalence between mass and energy, we use these two words interchangeably. 

As will be explained in Sec. 2.6, the masses calculated based on the brane geometry, are not 
really energies on the brane, from braneworld point of view. Actually they play no direct roles in 
braneworld, therefore we put “mass” in quotes in the title of this section. The purpose of studying 
these quantities is to compare the braneworld with GR to examine whether there are observational 
differences. In previous sections, the dimension of the spacetime is arbitrary. In this section, we 
only consider 3+1 dimensional spacetime (the brane). 

We will first introduce ADM mass and Hawking mass in Sec. 2.5.1 and Sec. 2.5.2, then we will 
derive the ADM mass in the conformally flat space (2.46). 

2.5.1 ADM Mass 

The ADM mass is introduced in standard texts [29, 32, 36]. One way to obtain it is to use 
(2.34) by setting the lapse a = 1 and the shift = 0 [36]. Note the ADM mass is only defined 
when the background is asymptotically flat. Applying (2.34) to a asymptotically fiat spacetime, 
the ADM mass reduces to [36] 

M A dm = 1 „ * lim < £ l&qij - % ( q kl qki )] s\ (2.55) 

167TCt4 Q—>oo J s Q 
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where qij is the spatial metric, and i,j = 1,2,3 is the spatial coordinate index. q,j is the (fiat) 
background metric, and is the covariant derivative associated with q\j. s 1 is the unit normal 
vector of Sf pointing outwards, which was denoted as Q a in (2.34). 

Note, only the spatial metric is needed in the definition. 

For the spherical symmetric case, the metric can always be rewritten as 


d / 2 = g rr dr 2 + r 2 (d 9 2 + sin 2 9d(/) 2 ) , 


(2.56) 


in spherical coordinate (r,0,<j>), where r is the areal radius. In this situation the ADM mass can 
be derived to be 

Madm = lim Afadm, (2-57) 

r—>oo 


where 


Afadm = 


r {g r r - 1 ) 


(2.58) 


2Cr4 \fgrr 

which is defined as the integrand in (2.55). In the special case of Schwarzschild metric where 
g rr = (1 — 2 mC? 4 /r) _1 , M a( j m = m/y/l — 2mGi/r , which has r dependence. Within the horizon, 
2mGi/r > 1, therefore this quantity is not well-defined within the horizon. 


2.5.2 Hawking Mass 

The Hawking mass is defined on two dimensional surfaces S with spherical topology, and it 
characterizes the mass in the space enclosed by the surface. The Hawking mass is described in terms 
of spin-coefficient formalism, for which please refer to [83, 85, 86 ]. I 11 terms of spin-coefficients, the 
Hawking mass is defined as 


M h = 


A 


I 67 tG 2 


1 + 2 nf s PP ' dS 


where A is the area of S. The p and p' are two of the spin-coefficients 


P = 


P = 


1 m af) X7 a (np + sp) = ^ e +> 


2 V 2 

2 V 2. 7 




V Q (np - sp) = 


2^2 


e_, 


(2.59) 


(2.60a) 

(2.60b) 
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where 0± are defined in (2.13). Substituting (2.60a) and (2.60b) back into (2.59), we get 


M h = 


A 


1 + 77^ 4> ©+9-d5 ) . 


167tG 2 V 167T./5 


(2.61) 


At an apparent horizon 0 + = 0, therefore the following mass-area relation holds at the apparent 
horizon 

/ A (2.62) 


Mh = 


167tG? 


In the case of spherical symmetric metric (2.56), the Hawking mass is calculated as 


Mh = 2 CU ^ ~~ ^ 9rr ' > ^ 


(2.63) 


where we have used 


0+ = = - V (fli 


which are directly calculated by applying the definitions (2.13) to metric (2.56). However, the 
calculation of 0’s depends on the time components of the metric too. This result is only valid in 
case the configuration is time symmetric (where the extrinsic curvature of t w constant hypersurface 
is zero). This motivates us to give (2.63) a new notation Mh, which agrees with Mu in the time 
symmetric case. 

In the special case of Schwarzschild metric where g rr = (1 — 2mG^/r)~ 1 , Mu = m and has no 
r dependence. 


2.5.3 The ADM “Mass” of the Conformally Flat Space 

Let us now apply the ADM mass or Hawking mass for the conformally flat space (2.46). First 
let us rewrite it in terms of the areal radius r, which is defined as the radius associated with the 
area of the r = constant surface. Therefore 


r 


= re 


A 


(2.64) 


Rewrite (2.46) in terms of fi, we have 

d^brane = (1 + rdA/dr)~ 2 dr 2 + r 2 (d# 2 + sin 2 6 d(j> 2 ) , (2.65) 
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which is the spatial metric on the brane. Substitute in the asymptotic behaviour of A at z = 

l,r> 1: A « ai/rz = aq/r, we have 

d^brane = (1 - aq/r) 2 dfi 2 + r 2 (d6' 2 + sin 2 6 &<j> 2 ) . 


Substitute this into (2.58) or (2.63), we obtain the brane ADM mass as 


AfbraneADM = CC1/G4. 

(2.66) 

Comparing this with eq. (2.54), we obtain a somewhat surprising result 


-^total — (^4/^5) ' -^braneADM- 

(2.67) 


2.6 A Quest for Brane Energy 


There is energy exchange between the bulk and the brane. To quantitatively describe the 
energy exchange, we need to define the energy in the bulk, and the energy on the brane. Basically 
only the energy in the bulk, or the energy on the brane, is sufficient since the other can be defined 
by subtraction from the total energy. 

To serve as the energy of the brane, an expression should satisfy: 

(1) its value is a part of the total energy. This rules out the ADM mass defined based on the 
brane geometry since it is equal to the total energy for a class of space configurations. Please 
refer to eq. (2.67). 

(2) it is not conserved during the evolutions, because of the energy exchange between the brane 
and the bulk. This requirement rules out any quasi-local definition on the brane, such as the 
ADM mass and Hawking mass evaluated from brane geometry. This is because, the quasi-local 
energies are defined as the limit at spatial infinities utilizing the time translational symmetry 
at spatial infinities, therefore conserved. 

The definition of an energy, especially that for the brane, is subtle and this section is to request 
the study of energy on the brane, rather than providing a solution. The reader can skip this 
section from here since the following is an attempt (that has conceptual issue) rather than a result, 
although the attempt turns out to have surprisingly good features exhibiting in the simulations we 
carried out. 
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We would like to search the definition from the Hamiltonian-Jacobi analysis described above. 
Since the energy of a certain region is defined as an integration over the boundary enclosing the 
region, if this concept could be generalized (which actually can not), we can tentatively define 
the energy of the bulk as the integration of (2.34) over the closed surface of Fig. 2.2(b). i.e. we 
enclose it along the brane. The difference, which is the integration of (2.34) along the brane, can 
be accordingly defined as the energy on the brane. 

However, without going into the details, we know this tentative definition is problematic for the 
following reasons. There is a sharp corner in the surface, where the extrinsic curvature is infinite 
and the integration of the extrinsic curvature over the corner can be indefinite. In the background 
spacetime, integrating along the brane makes sense physically. There is a requirement, however, 
that the surface is embedded into the background so the intrinsic metric of the surface stays the 
same, which can not be met for any non-trivial brane. Or alternatively, we keep the “embedding” 
requirement, but then we have to give up the “integrating along the brane” in background spacetime 
(which also means we give up the freedom fixing condition (2.39).). Also, since this is a definition 
regarding an integration on the brane, which will fail when there is a physical singularity on the 
brane. Therefore this definition can not be treated seriously. It is only a (very) rough description. 
We hope the resulting quantity changes monotonically with the amount of the energy exchanged 
between the brane and the bulk. 

In the simulation (presented in Chap. 5), we will choose the brane to be the surface in the 
background spacetime, i.e. we give up the requirement that the intrinsic metric of the surfaces 
in physical spacetime and background spacetime are the same. (2.42) and (2.43) are the two 
requirements for embedding, and we have to keep one of these to make the definition unambiguous. 
Here we choose to keep (2.42) because it is simpler and also because the areal radius still makes 
sense when there is a physical singularity (where condition (2.43) fails). 

Considering requirement (2.42), it is natural to express the energy in terms of spherical coordi¬ 
nate (r, 9, <j>) with f as the areal radius, where the physical metric is d l 2 = gffd^+r 2 ( dd 2 + sin 2 9dg ) 2 ). 
Since this definition only includes the integration on the brane, we may rewrite it as 

£bra„e = J a ^ k ) , (2.68) 

= dm f dr [ d# f d (j) v^det[q] a (( d ~ 2 ^k — . (2.69) 

Q-J-oo 87 tG 5 J o Jo Jo ' ' 

where f = f is the areal radius and f is the areal radius in the background space, q is the spatial 
metric of the brane (in physical spacetime) and we have y det [< 7 ] = r 2 sin 6 ^/gff. B stands for the 
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spatial part of brane manifold (in physical spacetime). Please note it is G 5 (rather than G 4 ) in 
the expression, and a is the lapse in 5D restricted on the brane, rather than the lapse function 
calculated from brane geometry. ( d ~ 2 '>k is the extrinsic curvature of the 3-brane embedded in 4D 
spatial bulk. Note a crucial feature of this energy is that, it is an integration over B, and its 
meaning is the energy over B. i.e. it might be possible to think of energy density, which is actually 
(spatial) coordinate independent on the brane. Noticing the second requirement of the energy on 
the brane (that the definition should not be semi-local), this definition might have captured a key 
feature. 
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Chapter 3 


Axisymmetric Spacetime with 
Non-Flat Background 


The cost of numerical computation increases dramatically with the number of dimensions. This 
is usually called the curse of dimensionality. To date, it is impractical to directly perform numerical 
calculations in 4 + 1 dimension at a reasonable resolution, or study problems with the demand of 
high resolution (such as critical phenomena) in more than (effectively) 2 + 1 dimension. Therefore, 
one of the first steps to study a physical system in higher dimensional spacetime, is to consider 
the system with some kind of symmetry, such as spherical symmetry and axisymmetry, to reduce 
the effective dimension. 

However, in many situations, singularities, instabilities, unavoidable noises or other irregulari¬ 
ties, arise in the numerical calculations performed under the coordinates adapted to the symmetry, 
either at the origin or at the axis of the symmetry ([50] and references therein). On the other hand, 
this kind of issue does not occur in the simulations of the same system performed under Cartesian 
coordinates. This kind of behavior is called the regularity problem in non-Cartesian coordinates. 
To make the discussion and presentation smoother, we will only mention axisymmetry below, 
which also applies to spherical symmetry. 

Where does the irregularity come from? Since the irregularity appears in cylindrical coordinates 
or spherical coordinates rather than Cartesian coordinates, it is widely believed that the irregularity 
comes from coordinate system choice ([50] and references therein). The terms involving l/r n (where 
r symbolically stands for the radius coordinate in spherical coordinates or cylindrical coordinates) 
often appears in the equations, and it is widely believed that, with or without machine precision 
playing a role, these terms are responsible for the irregularity [50]. 

What we do differently in this thesis is to make distinction between the coordinate system , and 
the fundamental variables —the unknown functions to be solved for—used in simulation (which are 
usually the components of tensors and pseudo-tensors), and reveal that the fundamental variables 
(rather than the coordinate system, or the operators in the coordinate system) are responsible 
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for irregularities. In particular, neither l/r ra -terms nor machine precision, is relevant to regularity 
problem. Therefore the regular results obtained from simulations performed in Cartesian coordi¬ 
nates, are actually due to the fact that the tensor (and pseudo-tensor) components are regular in 
Cartesian coordinates, rather than the regularity of Cartesian coordinates itself. As one will see in 
this chapter, actually it does not make sense to talk about the regularity of a coordinate system. 
The key to avoid irregularity issue in any coordinate system, is to construct tensors (and pseudo¬ 
tensors) in terms of regular variables that are compatible with the finite difference approximation 
(FDA) scheme at the vicinity of the axis (or the origin). There are many ways to express tensors 
(and pseudo-tensors) in terms of regular variables. To embody the construction, we present a 
general method to construct regular variables out of Cartesian components (which play the role 
of regular variables). The method can, for example, enable the generalized harmonic formalism 
to be used in non-Cartesian coordinates. Then we analyze why certain other formalisms in the 
literature can avoid regularity issue as well. 

Utilizing the knowledge obtained in studying regularity problem, the evolution schemes such 
as the generalized harmonic formalism and the BSSN formalism of GR can be rewritten under 
general coordinates which overcome the irregularity issue. On the other hand, there is another 
problem associated with the braneworld: the asymptotic spacetime background of braneworld is 
not flat, while the existing knowledge in the literature are only for asymptotically flat spacetimes. 
It is not clear how to setup the source functions (of the GH formalism) in the braneworld. The 
second part of this chapter is devoted to solving the non-flat background problem. 

3.1 Regularity Problem and Our Conjecture 

Very often (but not always), the numerical calculations performed in cylindrical coordinates 
adapted to axisymmetry, generate irregularities in the vicinity of the symmetry axis, but the same 
calculations performed in Cartesian coordinates do not generate irregularities. This phenomena 
is called the regularity problem. The irregularity, can appear as: (i) singularity in certain funda¬ 
mental variables used in numerical calculations in which situation the code would crash and the 
fundamental variables would not converge, or (ii) non-smoothness of certain fundamental variables, 
or (iii) smooth fundamental variables which do not converge at the expected order — i.e. they can 
not pass the independent residual tests. 

In this section, we will first introduce the existing analysis of regularity problem, then identify 
some of the existing methods that can overcome regularity problem, then we will declare a conjec- 
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ture regarding the true source of the problem. We back up the conjecture by a detailed analysis 
of the deviation between the numerical results and the analytical results, which reveals the second 
key element to yield regular results. 

3.1.1 The Existing Analysis and The Existing Solutions 

Let us take a specific example as the carrier of the existing analysis: the wave equation of a 
scalar field $ in 3+1 (flat) spacetime with axisymmetry. ( x,y,z ) are the Cartesian coordinates, 
and (p, cj), z ) are the cylindrical coordinates. The equation of motion is 

(V 2 - d tt ) $ = (d xx + d yy + d zz - dtt) $ = [d PP + ~^dp + 8 ZZ - dtt'j $, (3.1) 

where we have applied the fact that the space is axisymmetric so that the (^-derivatives are zeros. 
It is “obvious” that the term (1 /p)d p Q was singular and it was widely believed ([50] and references 
therein) that terms like this were responsible for regularity problems. 

Based on the analysis, the crucial step to cure irregularity was to modify the differential op¬ 
erators so that terms like 1 / r n did not appears. There are quite a few widely-used techniques to 
cure the regularity problem. Among these techniques, a class of methods called Cartoon meth¬ 
ods [50, 61], utilize the observation that “there is no regularity issue in Cartesian coordinates” to 
modify the operators in a way so that the operators are effectively Cartesian. 



Figure 3.1: The demonstration of the cartoon method. The system is axisymmetric with respect 
to 2 axis. For stencils which involve three grid points, only three slices are needed: the y = 0 slice, 
the one above it and the one below it. Note, this diagram is from [50]. 


The Original Cartoon Method 

The original cartoon method [50] was proposed by Alcubirre et al in 1999. Here I will again use 
a scalar field $ to demonstrate the basic ideas. The physical space is three dimensional space, with 
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(.x , y , z) being the Cartesian coordinates, and the system is axisymmetric with respect to z axis. 
Analytically, only one slice, say the y = 0 slice, is needed to include all the information because of 
the axisymmetry. Numerically, the cartoon method is trying to use only y = 0 slice as well. The 
method can be implemented via the following four steps. 

(1) Write the code using equations expressed in Cartesian coordinates. 

(2) Use only the y = 0 slice. Of course, this is impossible in Cartesian coordinates, since the 
discretization of derivatives with respect to y needs more than one slices in the y direction. To 
be specific, let us assume the equations consist of second order derivatives at most, and the 
discretization is of the second order accuracy. In this situation, the stencil only needs three 
grid points to do the finite difference approximation for the differential operators (Fig. 3.1). 
However, the function values on the upper and the lower slices are not known. 

(3) Obtain the function values on the upper and the lower slices by numerical interpolation con¬ 
structed from the values within y = 0 slice, utilizing the axisymmetry. The process is shown 
in Fig. 3.1. 

(4) Replace the function values on the upper and the lower slices, by the interpolated values 
obtained in step (3), and then substitute into the discretization. 

From these four steps, operators such as l/r" are avoided, and only Cartesian operators are directly 
used. 

The Lie Derivative Cartoon Method 

In the third step of the original cartoon method, in order to get derivatives with respect to 
y, numerical interpolation is used. The Lie derivative cartoon method [61] improves this, by 
analytically replacing these derivatives by those within the y = 0 slice. For the particular example, 
using the fact xd y — yd x being the Killing vector, we have 

0 „* = -<>:*■ 
y X 

Taking derivatives with respect to y on both sides, we have 

dy V $ = -d x § + -d xy $. 
x x 
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The laplacian V 2< h can then be evaluated at the y = 0 plane as 


V 2 ^ = (d XX + dyy + d ZZ ) 4> 



H— d x + —d xy 
x x 



d> 


y =o 



+ —d x + d zz \ 


T. (3.2) 


The derivatives in y direction are replaced by the derivatives within the y = 0 slice, on the analytical 
level (rather than numerical level in the original cartoon method). 


Background Removal Method 

The cartoon method is a general approach to solve regularity problem, which can be used in 
any formalism of GR. On the other hand, Brown [97] and Gourgoulhon [36] developed a method to 
rewrite the BSSN and generalized harmonic formalisms, by a background removal method, which 
will be described in detail in the following sections. The same method was also applied in [21] 
to solve static problems using Ricci-DeTurck flow methods. The simulation using this rewritten 
BSSN formalism in spherical coordinates “turned out to be” regular [56]. However, the authors 
did not analyze why this method is regular. Sorkin-Choptuik [54, 55] used a different method 
to remove the background effect, which again “turned out to be” regular, without analyzing the 
reason. 


3.1.2 The Conjecture: Variables Rather Than Coordinates 

Various other attempts on the solution to regularity problem had been proposed (see [50] and 
references therein) and they were not very successful until the cartoon method appeared. The 
cartoon method is a general method that solves the regularity problem, and there exist other 
methods which overcome the regularity problem as well[36, 71, 97]. Here we will take a closer look 
and provide a deeper understanding of the topic. The understanding can serve as a criteria and 
a guideline to develop regularized formalisms, and can potentially completely solve the regularity 
problem associated with coordinates. 

We start from the facts rather than speculations: 

(1) The Lie derivative cartoon method works well [61-63]. 

(2) 1/x still appears in eq. (3.2) which was produced based on Lie derivative cartoon method. 

(3) Eq. (3.2) is the same as (3.1) (equation in cylindrical coordinates), if x is identified with p. 

(4) Certain background removal methods (performed in non-Cartesian coordinates ) [36, 54, 55, 97] 
are free from regularity issues. 
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(5) There exist many successful simulations with l/r n terms (many references, e.g: [59, 61-63, 71]). 

(6) All in all, a multiplication by 1 /r n terms does not change the numerical feature. (p ■ d pp + d p ) $ 

is not problematic, therefore ( d pp + d p /p ) should be as safe. 

These facts lead me to ask whether the spherical/cylindrical coordinates were the sources of 
the problem, and ask the question: does the regularity issue come from the terms (1 /r n ) in 
the operators associated with the coordinate system, or somewhere else such as the fundamental 
variables being used in numerical simulation? 

In fact, if we take a closer look at the difference between the past simulations in cylindrical 
coordinates and those in Cartesian coordinates, the situation is either to use cylindrical components 
(of tensors and pseudo-tensors) under cylindrical coordinate, or to use Cartesian components 
under Cartesian coordinates, i.e. whenever one switches from cylindrical coordinates to Cartesian 
coordinates, he also “naturally” switches the fundamental variables from cylindrical components 
to Cartesian components, i.e. he has performs two changes: change of coordinates, and change of 
fundamental variables. It is not clear whether the regular results obtained from simulations under 
Cartesian coordinates, are due to the change of coordinates, or due to the change of fundamental 
variables. All the previous studies, including those yielding regular results, did not make the 
distinction between the effects from coordinates and the effects from fundamental varibles. Here we 
will make the distinction and make the conjecture: The regularity issue come from the fundamental 
variables used in simulation, rather than the operators associated with coordinate systems. 

According to the conjecture, the fact that simulations in Cartesian coordinates with Cartesian 
components are regular, is because the Cartesian components (of tensors and pseudo-tensors) are 
regular, rather than “Cartesian coordinates are regular”, or “operators in Cartesian coordinates 
are regular” 9 . According to the conjecture, the key to avoid regularity issue, is to express various 
tensors and pseudo-tensors in the equations in terms of regular functions 10 . In the following, we 
will first construct regular functions out of Cartesian components, then construct an example as a 
direct test to our conjecture. 

9 According to the conjecture, the fact that the cartoon methods produce regular simulations, is because the 
Cartesian components (of tensors and pseudo-tensors) are used as fundamental variables, rather than that the 
operators are made effectively Cartesian. 

10 According to the conjecture, eq. (3.1) should be free from regularity problem. The simulation confirms this 
claim. 
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3.1.3 Cartesian Components Method 

It is fairly easy to express a tensor in terms of regular components and there are multiple ways 
to achieve this. As a specific example, here we adopt a “canonical” approach: since the Cartesian 
components are regular, one way to express tensors in terms of regular components, is to rewrite 
tensors in cylindrical coordinates in terms of the (regular) Cartesian components, using the basic 
tensor transformation relation. Since there is a geometrical relation between Cartesian coordinates 
and cylindrical coordinates, we will do the coordinate transformation at certain region of the space, 
such that the radius coordinate p (or r) coincides with one of the Cartesian coordinates (refer 
to Fact 3 in Lie derivative cartoon method above), and the functional forms in the cylindrical 
coordinates, are the same as that in the Cartesian coordinates. Let us call this procedure as 
Cartesian component method. 

To clarify the concept and demonstrate the usage, let us take metric function in 3 +1 dimension 
as an example. We use indices “cylin” and “Cart” to indicate cylindrical and Cartesian, respec¬ 
tively. The procedure to express tensors in terms of their Cartesian components in cylindrical 
coordinates, comprises the following four steps. 


(1) write down the most general form of the metric (according to axisymmetry) in cylindrical 
coordinates (t, p , <f>, z) 


(cylin) _ 


( .(cylin) 
Utt 

(cylin) 

Utp 

0 

(cylin) \ 
3tz 

(cylin) 

<h P 

(cylin) 

ypp 

0 

(cylin) 

ypz 

0 

0 

(cylin) 

Vu 

0 

(cylin) 

\ i Hz 

(cylin) 

ypz 

0 

(cylin) 
!Jzz / 


(2) apply coordinate transformation 


(„) _ dx a dx? (x) 
9ttv dy 1 * dy v 9oc ^ ’ 


(3.3) 


(3.4) 


to transform this metric into its Cartesian coordinates (t, x , y, z) at the location (y = 0, x > 0), 
which is f) = 0 half plane in cylindrical coordinates, where the positive half of the x axis 
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coincides with p. Now the metric in Cartesian coordinates at (y = 0, x > 0) is 


1 mt Vtx 

0 

Viz ^ 



Pl.x Pxx 

0 

Vxz 



0 0 

Vyy 

0 



y Vtz Vxz 

0 

Vzz ^ 



( (cylin) 

gtt 

(cylin) 

9tp 


0 

(cylin) \ 
9tz 

(cylin) 

fJtp 

(cylin) 

9pp 


0 

(cylin) 

gpz 

0 

0 

4f"7* 2 

0 

(cylin) 

y 9tz 

(cylin) 

9pz 


0 

.(cylin) i 
gzz y 


where the first matrix is defined by 


(3.5) 


(3.6) 


y V (t, x,z)= £^ art) (t,x,y,z)\ 

ly=0,;r>0 

which is merely the restriction of Cartesian components onto (y = 0, x > 0) half plane. The 
second matrix is the calculation result of (3.3) via coordinate transformation relation (3.4). 


(3) compare the components in (3.5) and (3.6), and rewrite the tensor in terms of Cartesian 
components. Take the tt component as an example: the above relation tells us that x, z) = 

(i> P , z ) f° r all p = x > 0, therefore ryt. and have the same value at every (t, p , z), 

which means they have the same function form in terms of x (and p). Since x , y, z ) is 

regular, its restriction to (y = 0, x > 0), x, z), must be regular, which means (t , p, z) 

is regular in cylindrical coordinates. Actually all functions are of the same form as the Cartesian 
components, except for <fxj> component, for which we have g^ hn ' > /x 2 = r\ vy (and remember, 
x = p), and this relation suggests to rewrite the (f>4> component as p 2 y yy . 

(4) Finally we assemble components together to express the metric in cylindrical coordinates in 
terms of Cartesian components. The metric in cylindrical coordinates can be written as 


g(cylin) _ 


Vtt 

Vtp 

0 

Vtz 

1 

Vtp 

Vpp 

0 

Vpz 


0 

0 

p 2 1U 

0 


Vtz 

Vpz 

0 

Vzz 

J 


(3.7) 
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which is the same as 


Vtt 

Vtx 

0 

Viz 

Vtx 

Vxx 

0 

Vxz 

0 

0 

P 2 Vvv 

0 

Viz 

Vxz 

0 

Vzz 


with the subindices renamed. In the future we will denote this process roughly as 


(3.8) 


V^{t,p,z) = g 


(Cart) 


(t,p,0,z). 

y=ij=(p,x=p 


To this point the procedure of “expressing (metric) tensors in terms of their Cartesian compo¬ 
nents” is complete. Afterwards, Einstein’s equations (and other tensorial equations) are expressed 
in cylindrical coordinates to perform numerical calculations, which are regular. 

Let us study how this method is related to Lie derivative cartoon method. The residuals of 
Einstein’s equations are = G The following two equations are obtained by direct 
calculation, via Lie derivative cartoon method and our Cartesian components method, respectively. 


( ^ DC 

34 DC 

0 

3& DC N 

34 DC 

»LDC 

, '-xx 

0 

^? C 

0 

0 

»LDC 

Jl yy 

0 

V 8& DC 

»LDC 

Jv xz 

0 

^?° , 


^ t C 


0 

«£ c ) 


( 

3& DC 0 

5R)f c 



0 



54 DC 

^ DC o 

^? C 

0 

0 

»cc 

0 


0 

0 p 2 Ryy C 

0 



0 

) 


v *fc DC 

^DG 0 

5RL? C 


(3.9) 


(3.10) 


where LDC stands for Lie derivative cartoon method, and CC stands for our Cartesian components 
method. LDC is expressed under Cartesian coordinates, while CC is expressed under cylindrical 
coordinates. The result means, with identifying x = p, the two residuals are the same, except 
for a multiplication of p 2 in the <j>(j> components in cylindrical coordinates. However, as explained 
in fact (6), or Sec. 1.6, a multiplication of a smooth, non-zero function onto a residual equation, 
does not change the numerical properties of the numerical calculation. Therefore, our Cartesian 
component method is the same as Lie derivative cartoon method via a different approach under 
different philosophy. 

This agreement is not surprising. We can actually “prove” it as the following: both the Lie 
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derivative cartoon method and our Cartesian components method, are expressing the same tensor 
(the Einstein tensor) under the same coordinates (with identifying x with p), at the same physical 
location (y = 0 « (j>, x = p > 0), using the same variables (Cartesian components), therefore the 
results from the two methods must be the same. 

In the end, we reiterate that there are multiple ways to express tensors in terms of their regular 
components (such as the example with a superficially singular metric in Sec. 3.4, or background 
removal method in generalized harmonic formalism and BSSN formalism to be introduced below). 
Here expressing the metric (and other tensors or pseudo-tensors) in terms of their Cartesian com¬ 
ponents using the transformation relation of tensors and pseudo-tensors, is just a specific example 
to obtain regular components. 

Also, the above example can actually be improved: using local flatness [52] at p = 0, the above 
ri^ can be replaced by y pp + pW with W\ p= o = 0 (or replaced by r\ pp + p 2 W with W >p | p= o = 0). 
The local flatness condition can also be obtained by the procedure in Sec. 3.1.4. 


3.1.4 Results for General Symmetric Tensors 


In GR, the governing equations G pl/ = kdT pu are symmetric tensors of (0, 2) type (i.e. tensors 
with two “downstairs” indices). As long as we know how to deal with symmetric tensors of (0, 2) 
type, we know how to deal with Einstein’s equations — we know how to deal with GR. 

Let us consider a general symmetric tensor of (0, 2) type which in general has the following 
expression in Cartesian coordinates (t, x, y, z ) 


^i(Cart) 


T tt 

Ttx 

Tty 

T tz } 

Ttx 

T xx 

Txy 

T 

X xz 

Tty 

T 

± xy 

T 

± yy 

Tyz 

T u 

T 

X xz 

T 

± yz 

T zz ) 


(3.11) 


Now we will express T in cylindrical coordinates. Since the final expression should be independent 
of 0, one can express T at any value of (j>. e.g. let us choose <f> = 0. At this “location” we have 
x = p, y = 0. 11 

The axisymmetry is expressed as 


£{-yd x +xdy)T — 0 . 


(3.12) 


11 Again, we emphasize that the expression in cylindrical coordinates does not depend on the value of <f>, and here 
we just take advantage of this fact and do the calculation at (f> = y = 0. 
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Opening up this expression, we obtain the following relations at y = 0 


Tty — T xy — T yz — 0 , 

Ttt.y = T tx ,y = T tZt y = T XZ y = T XXt y = T zz< y = Tyy^y = 0, 

T tx = x ■ U, T xz = X ■ V , Tyy - T xx = X -W, where W| x=0 = 0, 


(3.13) 

(3-14) 

(3.15) 


where (U, V. IV) are regular expressions in terms of the T pu , which specific forms are irrelevant at 
this point. 

Performing a coordinate transformation from Cartesian coordinates (t, x , y , z) to cylindrical 
coordinates (t,p,<p,z) and using (3.13) (and take <j> = 0), we obtain 


'J’(cylin) 


Ttt 

T tp 

0 

T tz 

Ttp 

T pp 

0 

T PZ 

0 

0 

P ■ Tfifr 

0 

Ttz 

T PZ 

0 

Tzz 


(3.16) 


where r MJ/ (t, p, z) = T fll/ \ y= ^ =Q x=p (t, p, 0, z), which are guaranteed regular if the Cartesian com¬ 
ponents are regular. 

The condition (3.15) now reads 


Ttp ~ P U, 

(3.17) 

Tpz= P- V, 

(3.18) 

t^ = T pp + p ■ W such that W \ p= o = 0. 

(3.19) 


The first two can be alternatively expressed as Tt p \ p= o = r pz \ p= q = 0 12 , which are parity conditions, 
and the last condition is called local flatness. To complete the parity conditions at p = 0, similar 
to the derivation of (3.14) at y = 0, now opening up (3.12) at x = 0, and then renaming the indices 
and notations, we get 


T tt.p|p=o ~ T tz,p\p=o - T pp’p\ P =o - T <t><l>’p lp=o ~ Tzz A p =o - (3.20) 

i.e. the local flatness condition and the parity conditions about the metric tensor obtained in 
Sec. 3.1.3, actually apply to any symmetric tensors of ( 0 , 2 ) type. 

12 For functions that can be expressed as Taylor expansion f = fir 1 , the condition f\ r =o = 0 implies /o = 0. 

Therefore / = YaL\ fi r * = r ' T,i= o fi+ i r ’ = r -V. 
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3.1.5 Test of the Conjecture 

Here we study a testing problem: massless scalar field evolution in 5D. By the same procedure 
as what is used in 4D, the metric in cylindrical coordinates ( t,r,Q,<j>,z) is 


& 


Vtt 

Vtr 

0 

0 

\ Vtz 


Vtr 0 0 Vtz ^ 

Vrr 0 0 Vrz 

0 veer 2 0 0 

0 0 Veer 2 sin 2 9 0 

Vrz 0 0 Vzz ) 


(3.21) 


Einstein’s equations in terms of (3.21) are equivalent to that from Lie derivative cartoon method, 
and the simulation is, without surprise, regular. 

Now, to directly test our conjecture, we also perform the simulation using the following metric 
representation 


^singular 


f Vtt Vtr 0 

Vtr C ^ 0 

0 0 Veer 2 

0 0 0 

^ Vtz Vrz 0 


0 Vtz ^ 

0 Vrz 

0 0 

Veer 2 sin 2 9 0 

0 Vzz 


{3.22) 


i.e. we purposely introduce a “singular” term £• (r+ l)/r (or, we have defined £ = yxy • Vrr): which 
would be troublesome from the conventional point of view. However, if our conjecture is correct, 
then the simulation in terms of (3.22) would be just as good as the simulation in terms of (3.21), 
since £ = • Vrr is regular. 

To test the claim, we performed the simulation in terms of both metric representations. It turns 
out that the results of the two simulations are the same (which means the superficially singular 
term £ • does not cause trouble). Therefore, this simulation supports our conjecture. 

For the details of the simulation, please refer to Sec. 3.4. 


3.1.6 Validity and the Extension of the Conjecture 

The conjecture is simple, and its usage and validity have been demonstrated though examples. 
However, the demonstration is not a proof. In this section we try to give a proof (in a physicist’s 
sense) and we will address the following three questions: 

(1) Why terms like l/r ra do not cause singularity problems? 
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(2) When taking only regular functions, are there still any other kind of regularity issues (other 
than the singularity issue)? And if yes, how to cure it? 

(3) What role does machine precision play in the regularity problem? 


To answer these questions, we study at what degree the finite difference approximation (FDA) 
can represent a derivative operator. It is clear to use an example. For a function /( x) with h as 
the spacing of the grids, let us consider the following FDA 


fi +1 fi— 1 

2 h 


= fix) + h 2 C 2 f" + h 4 C 4 f"'" + --- = Y, h 2 n C 2 n f {2n+1 \ 


n =0 


where C n are constants which specific values are irrelevant here, and /^ is the n-th derivative 
of /. The expression shows the FDA of a first order derivative, discretized at a finite difference 
scheme that is of the second order accuracy. Generally, the FDA (.4) of a d-th derivative operator, 
discretized at a finite difference scheme that is of o-th approximation order, is 


OO 

A f(d) = f(d ) + h a Ca f(a+d) + h 2a C2a f(2a+d) + . . . = ^ h na Cna f(na+d), 

71=0 


(3.23) 


where h is the grid spacing. 

In GR, when r is small, the leading orders of the Cartesian components of symmetric tensors 
of (0,2) type (denoted as /) are asymptotically either linear or quadratic in r (refer to Sec. 3.1.4). 
Therefore, as long as a + d > 3, we have f( na + d ) = 0 (n = 1,2,...) in (3.23), therefore 


A f(d) = f{d) 


(3.24) 


i.e. the FDA approximation with a + cl > 3, is exact when / is linear or quadratic in r. In general, 
when the FDA is exact in the vicinity of the symmetry axis, we say that the FDA is compatible 
with the asymptotic behaviour of fundamental variables. Therefore, as long as the equations are 
well-behaved at the continuous limit, they are also well-behaved at the discrete level. In particular, 
terms with multiplications of l/r n are well-behaved since the discretization results of the functions 
are exact, and multiplication operation is (almost) exact as well. 

The above equation explained why our solution to regularity issue works, and it can also predict 
the failure of FDA. From equation (3.23) one can see that if f( a + d ) ^ o, then /^ and h a C a f ( ' a+d 
are comparable at small r (because (a + d) — a = d), then the FDA is a bad approximation, i.e. the 
FDA has to be exact , otherwise the error {h a C a f^ a+d ' > ) is of the same order as the value (/W), 
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which invalidates the FDA. For example, when / = r 3 and a = 2, then the FDA for /' has trouble. 
Given the asymptotic behaviour of / being / ~ r 3 , one can find ways to cure, such as 

(1) Use higher order FDA. 

(2) f —» diff(/, r 3 ) x 3r 2 , then only discretize the “diff”. Here the notation “diff” is the partial 
derivative: diff(A,H) = dA/dB. 

(3) Rewrite / as / = r • F and let F be the fundamental variable instead. The criterion to rewrite 
/ in terms of F is that the asymptotic behaviour of F is compatible with the FDA in the 
vicinity of the symmetric axis. 

These method can be easily generalized to other asymptotic behaviours (for example, using the 
third method, a function with asymptotic behaviour as / ~ ?’ 5 / 2 can be rewritten as / = r 1 / 2 • F). 
The first two methods have been used in the literature and proved to be successful. Overall, our 
opinion is that, as long as the asymptotic behaviours are known, there is always a way to easily 
fix these issues. 13 

We end this section by analyzing machine precision and show that it is not related to the irreg¬ 
ularity issues. In the above analysis, we only analyzed truncation error (caused by the finiteness 
of h) with machine precision e (~ 10 -16 in double precision) being ignored and we concluded that 
both l/r n terms and non-(1 / r n ) terms are exact, as long as the FDA scheme is compatible with 
the asymptotic behaviours of the fundamental variables in the vicinity of r = 0. Now we restore e 
and analyze its effects. The concern regarding e is that the l/r n terms would amplify the effect of 
machine precision so that the errors in 1 /r n terms might be significantly larger than the errors in 
non-(l/r”) terms. In the following, we will show that the errors in l/r n terms are not larger than 
the errors in non-(l/r ra ) terms. 

For specificity, let us consider the FDA of 


f,rr + f/r 2 - 1/r 2 , 

13 However, in the literature there is yet another type of “cure” to regularity problem. Let us take / = r 3 using 
second order FDA as the example. The basic idea of the cure is to still discretize f as f —>■ (fi+i — fi-i)/2h, 
which is actually wrong according to our analysis above since the truncation error is comparable to the value itself. 
But it is easy to prove that, for any fixed r, the truncation error convergences to zero as the resolution increases. 
However, for any resolution, there is always a region in the vicinity of r = 0, where the error is still comparable to 
the value itself. In particular, the values at the first few grids in the vicinity of r = 0 will always be wrong. Since the 
mistake only happens in a finite region (which size shrinks as the resolution increases), the result in other regions is 
still correct, as long as the simulation is stable. Therefore the focus of this kind of “cure” is to make the simulation 
stable. Our opinion is that, these methods start from a mistake and try to do something to cover the mistake so 
that it is controllable within a shrinking region (as the resolution increases), and even if the mistake is controllable, 
the behaviour at r = 0 is always problematic (for example, the independent residual tests would not be passed at 
the point). Instead of searching for any “cure” of this kind, our suggestion is to develop methods based on our two 
criteria (see Sec. 3.1.7) on regularity, such as the three ways presented above. 
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where the asymptotic behaviour of f(r) in the vicinity of r = 0 is / ~ 1 + r 2 . The concern is that 
f rr ~ ~ 1, but f/r 2 would be of order 1/e when r is as small as ~ yTF, therefore the addition 

operation in f <rr + f/r 2 is not accurate. Symbolically, we denote this concern of non-accuracy as 

f,rr + f/r 2 ~ 1 + 1/e = e -1 • (1 + e), (3.25) 

whose accuracy is the same as 1 + e since multiplication operation does not lose accuracy. We 
will analyze what it takes to let this non-accuracy in f rr + f/r 2 show up. For the expression 
f,rr + f/r 2 — 1 /r 2 , the discretized equation at the second order FDA is (fi+i — 2/* + /,_i)//i 2 + 
(/— 1 )/r 2 . Representing / >rr as (fi+i — 2fi + fi-i) /h 2 , relies on the validity of Taylor expansions 
such as fi +1 = fi + hf t r + {h 2 /2)/ jrr + (/i 3 /6)/, rrr + .... It is then crucial to examine whether 
the higher orders in these expansions can be expressed accurately by floating point numbers. 
To let (fi+i — 2 fi + fi- 1 ) /h 2 be of second order accuracy, it is required that (/i 3 /6)/, rrr in the 
expansion must be representable by floating point number. Even if we give up the desire for second 
order accuracy, it is required at least (h 2 / 2) f^ rr is representable in the expansion by floating point 
number, i.e. it is required the value of (h 2 /2)f rr is not lost in /, + hf, r + (h 2 /2)f rr , which can 
be symbolically expressed as ~ 1 + \fe + e. Therefore in ( f t +i — 2/,; + /;_i) //i 2 , the operation is 
symbolically (1 + \fe + e — 2 + 1 — + e)/e. The accuracy level is determined by the largest 

and the smallest values involved in the addition (subtraction) operations, which are 1/e and e/e. 
i.e. the accuracy level of this FDA is 

e _1 • (1 + e). (3.26) 

i.e. in order to let the FDA (/j+i — 2 fi + fi-\) /h 2 represent / >rr at the lowest order accuracy, the 
e is required to be representable in the operation 1 + e. Comparing (3.26) with (3.25) one sees 
that, the accuracy level (due to machine precision) of non-(l/r n ) terms is the same as that of the 
l/r ra terms. If the operations involved l/r n terms (such as / >rr + f/r 2 ) lose their accuracy due to 
machine precision, then the FDAs of non-(l/r n ) terms can not represent the derivative operators 
any more (such as (fi+i — 2 fi + fi-i) /H 2 fails to represent f >rr at any grid point, rather than 
merely a few grids in the vicinity of r = 0). i.e. what it takes to fail f rr + f /r 2 , is the failure of 
FDA method. 

Or equivalently, provided the FDA is still valid (so that finite difference method can be used), 
then / rr + f/r 2 does not lose its accuracy (therefore machine precision can still be ignored), 
i.e. for the aspect of machine precision, the l/r n terms are not worse than the non-(l/r n ) terms. 
Therefore machine precision plays no role in regularity. The discussion above reveals that the 
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truncation error (due to the finiteness of h) is fundamentally different from the machine error. 
The former is a generic feature associated with the FDA method which has nothing to do with 
machines, while the latter is due to the limited ability of machines. One can improve the latter by 
using more accurate numbers such as quadruple-precision floating point numbers with the machine 
precision e ~ 10 -34 , while the truncation error can not be improved in this way. 

3.1.7 Summary of Regularity 

This subsection is the summary of the regularity discussion, which should be treated as a 
rephrase of our conjecture: the irregularity has nothing to do with coordinate system (and operators 
in the coordinate system) or machine precision. The simulation in cylindrical coordinates (or any 
other coordinates) is regular, provided the following two conditions are met: 

(1) The fundamental variables being used in numerical calculations, are regular functions; 

(2) The FDA scheme is compatible with the asymptotic behaviours of the fundamental variables 
at the symmetry axis (or the origin or any other special locations in the coordinate system). 


3.2 The Generalized Harmonic Formalism in 
Non-Cartesian Coordinates 

In the previous sections, we described how to apply the Cartesian components method to 
express the metric functions in terms of regular functions. In this section, we apply the same 
method to the generalized harmonic formalism (and to BSSN in appendix A). 


3.2.1 An Introduction to the Generalized Harmonic Formalism 


Now we briefly introduce the generalized harmonic formalism [61] of GR. The source functions 
II 11 are defined as 

H a = V p V p x a = -F“^, (3.27) 

and H /t = . In terms of H M , Einstein’s equations 


Rni/ — kd 


T _ 

^ d-2 


9tiuT 
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reduce to [61] 

-\a af) 9^ - g a \,,9„ w , a - H M + H^T\ V - r%r^ Q = k d (t^ - . (3.28) 

The generalized harmonic formalism [61] is to lift If, as fundamental variables , then the defining 
equations (3.27) are not defining equations any more. Instead, they are now constraints 

C a = H a - T“ v g l * v ~ 0, (3.29) 

where ~ means the equations are constraint relations, i.e. in the generalized harmonic formalism, 
the fundamental variables are the metric functions g^ v and the source functions , where If , are 
subject to the constraints (3.29). Since If, are now fundamental variables, the principal parts of 
(3.28) are now —^g a ^g fiu ,ai 3 , which are manifestly strongly hyperbolic. 

Using the Bianchi identities, it can be proved [61] that the constraints (3.29) will always be 
satisfied during an evolution, as long as they are satisfied at any instant of the evolution and the 
Hamiltonian constraint together with the momentum constraints are also satisfied at the instant. 
In numerical simulations, however, there are always modes violating the constraints, and the 
violation modes might grow with time (even exponentially), which will produce unphysical results. 
Gundlach et al. then suggested [58] to add constraint damping terms to the left hand side of 
eq. (3.28), which can be 

Zftu = k (n^Cv) - ^-^-g^nPCp] , (3.30) 

where k > 0, — 1 < p < 0. If the constraints are satisfied, the damping terms vanish. If small 
violation modes develop, the damping terms can damp out the modes [58] and drive the numerical 
results back to physical configurations. 

Since If , are now fundamental variables, their equations of motion need to be imposed. The 
constraints H a ~ — T“ v g 111 ' = show that the If, are related to coordinate gauge choices. 

There are a lot of freedom to impose these equations of motion. For example, If , can even be 
functions of coordinates and metric functions such as Hf = g^x 1 '. But generally it is required that 
If , do not include the derivatives of metric functions (which means the equations of motion such 
as H t , = gfj, u ^ a g ua are generally not recommended), so that the principal parts of (3.28) are not 
affected. We do not intend to extensively discuss the gauge choices in this section. Please refer to 
Sec. 5.3 for the details of some popular gauges being used in the literature. 
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To simplify future discussions, we define 


pa _ pa pv 
x x nvy 


3.2.2 The Generalized Harmonic Formalism in Cylindrical Coordinates 

Numerical evolutions using GH in cylindrical coordinates in terms of H are unsuccessful due 
to the fact that the // M 's are singular. Instead, our approach is again to use Cartesian components 
method: employ the coordinate transformation relation for Christoffel symbols 


(j,) r « _ d^dx v (x) dy K dy K d 2 xT 
a/3 dy a dyP dy7 + dx» dy a dyP ’ 


(3.31) 


to express U M ’s in terms of their Cartesian components. By the same procedure, we obtain 

ht + (1/p) (ntp/w) 
h P + (l/p) {Vpp/Vh) 

0 

hz + (1/p) (Vpz/vh) 

where h a (t,p,z ) = (Cart) R a (t,x = p,y = 0 ,z) with sub-index renamed, therefore are regular. 
Specifying the gauge amounts to choosing appropriate h a (rather than H a ). e.g. the harmonic 
gauge now reads h M = 0, which is indeed the harmonic gauge in Cartesian coordinates. From 
(3.32), one can see why the simulations with H M as fundamental variable fail—because H r is 
singular, which violates our first criterion for regularity. 

By the same procedure, the metric in spherical coordinates (t, r, 0, </>), with spherical symmetry, 
is 

rjtt Vtr 0 0 

r)tr rjrr 0 0 

0 0 peer 2 0 

0 0 0 peer 2 sin 2 6 
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Again, yee is replaced by r/ rr + rW , considering local flatness. The source functions are 




^ h t + (2/r) (vtr/vee) ^ 

H r 


h r + (2/r) (vrr/yee) 

H e 


cot 9 

V H <’ ) 


l 0 ) 


For 5D, the metric in cylindrical coordinates ( t,r,0,(f>,z ) is 


Cn 5 = 


\ Vt 


Vtt Vtr 
Vtr Vrr 
0 0 

0 0 

V rz 


0 

0 

Veer 2 

0 

0 


0 Vtz ^ 

0 Vrz 

0 0 

Veer 2 sin 2 6 0 

0 Vzz ) 


(3.34) 


(3.35) 


where vee is replaced by y rr + rW, considering local flatness. And the source functions are 


H t \ 


( h t + (2/r) (vtr/vee) ^ 

H r 


h r + (2/r) (vrr/vee) 

Hg 

= 

cot 9 

h 4> 


0 

Hz ) 


^ h. z + (2/r) (vrz/vee) J 


3.2.3 Background Removal in the Literature 

In the literature, there exist a class of methods for solving the regularity issue of the H^s 
through the use of the source functions with a background term subtracted, which can all be 
called background removal methods. 

In [54, 55], the authors used (for example in the axisymmetric case in 5D) 


H r = H r + -,H t = H t ,H z = H s , 


where 2/r is the value of H r in flat spacetime background, and the variables with a hat (') are 
the variables the authors used as fundamental variables. From the discussion above, we obtain 





/ y ge - rW 
V Vee 



= hr 


2 IT 
Vee ’ 


80 




3.2. The Generalized Harmonic Formalism in Non-Cartesian Coordinates 


which is related to h r by a regular term, therefore is regular. Similarly, their 

H t = ht + VlL 

r r]g e 

and H z = h z + - —, 
r rjgg 

are regular too (because r]t r and ii rz are asymptotically linear in r as discussed in section 3.1.4). 
On the other hand, there is another way to remove the background [36, 97]. Instead of using 
s, the authors developed the formalisms to use 

h* - -<v Kp-Kp)g af) 

as fundamental variables, where the “bar” means that the associated quantities are for the back¬ 
ground. is actually a vector (while // ;J is not). Although no simulations using this method have 
been reported 14 , this method should be able to generate regularized simulations if the background 
is chosen correctly. The reason is as follows: if the background is chosen to be flat, in Cartesian co¬ 
ordinates (where = 0) we have H^ = — = h M , where h M is the Cartesian component 

of the source function, which is what is used in our Cartesian components method. Since H^ is a 
vector, its transformation from Cartesian coordinates to cylindrical coordinates at (y = 0, x > 0) is 
simply tf t (cylin) = tf t (Cart) = h t , I^ cylin) = ijj; Cart) = h r , and ff| cylin) = tfi Cart) = h z , while other 
components are zeros. Therefore H^ = h M holds in cylindrical coordinates as well. i.e. when the 
background is flat, the fundamental variables used by the background removal method in [36, 97], 
are the same as the ones used in Cartesian components method. Therefore, by our conjecture, this 
method has no regularity issues, as long as the Cartesian components are regular and the FDA 
scheme is compatible with the asymptotic behaviour of the fundamental variables. 

This formalism can solve the regularity issue (by choosing the background properly), and the 
source function is now a vector. Furthermore, it can remove any background (especially non-flat 
backgrounds) 15 . We will develop it further in Sec. 3.3. 

14 We will use a generalization of this method in our simulation of the braneworld. 

1 r ‘ [ lowever. when the background is not flat, the regularity of the source functions is not guaranteed, in which 
case the background needs to be analyzed in a case-by-case basis. 
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3.3 The Generalized Harmonic Formalism in Non-Flat 
Backgrounds 

The spacetime background for RSII is 

ds 2 = 4 ( - d t 2 + dr 2 + r 2 (d 6 2 + sin 2 9d<p 2 ) + d z 2 ^j. (3.37) 

For this background, the h’s from eq. (3.36) are 

(h t ,h r ,h z ) = ( 0 , 0 ,--^. (3.38) 

Normally setting up the gauge choices is to let take specific values or satisfy specific conditions. 
For example, in the case for the harmonic gauge, the are zeros. However, the spacetime 
background is now non-trivial, in the sense that h z is not zero for the background, and it is not 
clear how to let gauges, such as harmonic gauge, make sense. In other words, it is not clear how 
to setup h^. 

Generally, in the current literature, the gauge choices for GH are all in Cartesian coordinates 
with a flat spacetime as the background. In the last section, we have developed the GH in non- 
Cartesian coordinates, where we have used the assumption that Cartesian components are perfect. 
For the braneworld, however, Cartesian components might not be sufficient in the sense that it 
is not clear how to easily employ the existing gauge choices used in the literature, i.e. actually 
Cartesian components are not perfect and we need to find a way to “get rid of” the background. 

3.3.1 Tensorial Source Functions 

The first way is to use the “background removal” methods mentioned above. For example, 
following [97], we employ the following functions to serve as the fundamental variables 

K-~ 9^ (r^/3 - T v afj ) g aP , (3.39) 

where the bar ( " ) stands for quantities and operations associated with the background. As 
discussed in section 3.2.3, when the background is taken to be flat, in (3.39) is the same as the 
used in “Cartesian components” methods. Yet, (3.39) can be applied to non-flat background. 
Another good feature is that h M is now a tensor (vector), which geometrically makes more sense. For 
the braneworld, if (3.37) is taken as the background in (3.39), we have h a = 0 for the background. 
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However, still this can not be easily linked with the existing work in the literature where the 
background is flat. i.e. it is not clear what gauge should be given to the h a . 

3.3.2 Conformal Transformation 

Noticing the background metric (3.37) is conformally flat, we can do a conformal transformation 
= ( z 2 /(?)g tlI ,, to obtain g^ u whose background is flat. Actually we can use a general conformal 
factor 'I' in the conformal transformation as g = 4 '~ 2 g^v The only requirements on this 4' are: 
(a) 4 ' = l/z when the metric is the background metric; and (b) 4< goes to t/z asymptotically at 
spatial infinities. Or even more generally, we define conformal transformation 

g^u = q gf_iv, (3.40) 

where q is a numerical factor that can be set as any value for convenience. A tilde (~) is used for 
the quantities and operations associated with g^ u . For example, the Christoffel symbols are 

r“ MI/ = 2 9 a/3 {9pv,n + <5W 3,v - 9nv,p) ■ (3-41) 

Repeating the derivations in [32] or [36] while keeping q and d general, we obtain 

R^ = R^FR^\ (3.42) 


where 


R$ = - ^ V. In 4/ - V Q V“ In 4> 

+ q2[d ~ 2) (V M In *V„ In 4- - V Q In *V“ In ) . (3.43) 

Substituting this into Einstein’s equations 

R» v ~ \g^uR = k d T^ v (3.44) 

yy R^ = k d (t^ - ^ ^ 2 ^^) ' (3.45) 

we have 

= k d (t^ - - flW, (3.46) 
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where T = g l ' ,J T IJ i J . The equations are Einstein’s equations with a modified right hand side. 

Depending on how the left hand side is rewritten, one can proceed with either one of the 
following two approaches. The first approach is to define the source function as 

H* = -f \fg aS , (3.47) 


which results in the following GH formalism with a conformal function 


- 9 a \j v) p, a - + Hp fV - f%f^ Q = k d (r M „ - - i?W. (3.48) 


For braneworld simulation, this formalism is sufficient since the conformal transformation let the 
background of g^ be flat which enables us to borrow the existing results regarding gauge spec¬ 
ification in the literature. Eq. (3.48) can be solved using (3.35) as the metric and (3.36) as the 
source function (every quantity needs to have a tilde though). 

Alternatively, R^ u can be rewritten using tensorial source functions. This approach does better 
in the following aspects: it is not limited to the case where the background of g is fiat, and the 
source function forms a tensor. The background spacetime of g is denoted as g /1U . We further 
define 


f^tOL _ pa ■pQ' 

[XV - [XV [XV1 

(3.49) 

Af /3 = C p aS g aS , 

(3.50) 

h ' 3 = -Af ' 3 = - (f * a4 - fM g aS . 

(3.51) 


Repeating the derivation in [36] (for the BSSN derivation though), while keeping d general and 
without requiring g = g, we obtain the final equation as 

- \g a ^cS7g~g^ - V a ~g^V u) ~g ap ~ V { Ji v) + KC a ^ - C a ^C\ v 
= k d (t^ - -R^-RiV- (3-52) 

3.3.3 The Implementation of Generalized Harmonic Formalism 

The discussions in this subsection regarding the implementation apply to all the GH formalisms 
(3.28), (3.48) and (3.52). Here we use (3.52) as a specific example. 

The (t, fi) components of equation (3.52) are the Hamiltonian constraint and momentum con- 
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straints. Eq. (3.51) are extra constraints due to the introduction of h. In this sense, performing 
the evolution using all the components in (3.52), is a full-evolution scheme with Hamiltonian and 
momentum constraints being satisfied. Constraints (3.51) are driven to be satisfied by adding the 
constraint violation damping terms [58] into (3.52). Therefore we give it the name full-evolution 
with source driving gauge. This is what has been adopted in the literature. However, it is not very 
clear how the coordinates condition is imposed via h M (or even H M ). 

Alternatively, we can adopt the same strategy as that in BSSN (therefore we give it the name 
BSSN-like method): the lapse and shift functions are specified in the “ordinary way” (such as max¬ 
imal slicing, 1 + log slicing, Gamma freezing, etc), while the Hamiltonian/momentum constraints 
are left un-evolved. Also, an evolution equation for h M will be derived, rather than using h M to drive 
lapse and shift. Using T /3 ^ a = i d a In g (and f ,3 /3a = | d a lnjjj), we obtain V Q < 7 a/3 = h 13 — g a/3 V a F, 
where F = 4 In ( g/g ), which is a scalar. Taking derivative with respect to t, we have 

= -V t 5 Q/3 V Q ^ - g a0 V t S7 a g^ + V t V^F, (3.53) 

within which there are terms which will be replaced by quantities with first-order and 

zeroth-order in Vt, via eq. (3.52). As for the Hamiltonian/momentum constraints, the damping 
terms for their violation modes could be constructed. Beyond these, h^ + AT M ~ 0 are still con¬ 
straints (which are similar to the constraints f* ~ ^( d-1 >f 1 , k — in BSSN formalism). 

Note, we have not implemented a code using BSSN-like methods yet. Performing simulations 
using this method is part of our future plans. 

The BSSN method is widely used in numerical relativity, and we have simply borrowed its 
spirit, in deriving the generalized harmonic formalism in non-Cartesian coordinates in non-flat 
background. For completeness, we also generalized the BSSN formalism to non-flat background 
with non-Cartesian coordinates. However, since we will not use the BSSN in simulating the 
braneworld, the development of the BSSN method is put in appendix A. 

3.4 Evolution of Massless Scalar Field under Cylindrical 
Coordinates 

In this section, we study a model problem: the massless scalar field collapse in 5D under 
axisymmetry. The reason for choosing the collapse in 5D under axisymmetry is that this system is 
close to our main project — the massless scalar field collapse in the braneworld, and is therefore an 
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instructive model problem. We will use this section to demonstrate the use of the various methods 
developed in this chapter. Namely: we will test the Cartesian components methods of generalized 
harmonic formalism. Then we will test our conjecture by introducing a seemingly singular term 
into the metric, while the fundamental variables are known to be regular. 

3.4.1 Initial Data 

At the initial instant, we assume the metric takes the following form 

ds 2 =e 2A (- At 2 + dr 2 + r 2 (d<9 2 + dcp 2 sin 2 9) + dz 2 ), (3.54) 

under cylindrical coordinates (f, r, 9, <f>, z), with r = 0 being the symmetry axis. The range of the 
coordinates are of course r > 0, 0 < # < 7T, 0 < 0 < 27r, —oo < z < oc. In this simulation we will 
keep the spacetime symmetric about z = 0 plane, therefore the range of z (in the code) is [0, oo). 

The initial data is obtained by solving the momentum constraints and the Hamiltonian con¬ 
straint. We take time symmetric initial data, therefore the momentum constraints are satisfied 
automatically, which leaves only the Hamiltonian constraint to be solved. The Hamiltonian con¬ 
straint reads 


A- + A** + + (A r ? + (A *) 2 + ^ ((A ,) 2 + (A* 2 )) = o, 


(3.55) 


where <f> is the massless scalar field. We choose its initial configuration to be 


<E>(0, r, z) = a ■ exp 



(3.56) 


which is symmetric about z = 0. Here we purposely choose the spherically symmetric initial 
configuration, to test the quality of our axisymmetric code by examining whether a spherically 
symmetric initial data will remain spherically symmetric. For the testing, we take a = 0.1, cr = 
0.25, r 0 = 1. 

The time derivatives of spatial components of the metric are taken as zeros since the initial data 
configuration is time symmetric. The other metric components are taken trivially gtr = gtz = 0, 
but the time derivatives of gtt,gtr,gtz are chosen such that ht = h r = h z = 0 (ref. eq. (3.58)) at 
the initial instant [62, 63]. 
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3.4.2 Evolution 


The evolution is performed using our Cartesian components method applied to the generalized 
harmonic formalism (eq. (3.28)), therefore the metric representation is (3.35) 


Vtt 

Vtr 

0 

0 

Vtz 

Vtr 

Vrr 

0 

0 

Vrz 

0 

0 

Veer 2 

0 

0 

0 

0 

0 

Veer 2 sin 2 6 

0 

Vtz 

Vrz 

0 

0 

Vzz 


(3.57) 


where vee is replaced by r/ rr + rW. considering local flatness condition at r = 0, as discussed in 
Sec. 3.1.3. The gauge sources are (3.36) 




( h t + {2/r)(vtr/vee) ^ 

H r 


h r + (2/r) (vrr/vee) 

H 0 

= 

cot 9 

H^ 


0 

Hz ) 


^ h z + (2/r) (vrz/vee) J 


During the evolution, the gauge is chosen to be ht = h r = h z = 0, which is equivalent to harmonic 
gauge in Cartesian coordinates. 

A constraint damping term [58] is added to the left hand side of eq. (3.28) 




1 +p 
2 



(3.59) 


where k> 0, —l<p<0. k = 6 and p = 0 were chosen during the simulation. 

The simulation is performed in compactified coordinates R = r/(r + 1) and Z = z/(z + 1). In 
this way the spatial infinities, r = oo and z = oo are mapped to R = 1 and Z ml. Since we focus 
on the behaviour in the region that is close to r = z = 0, where the resolution in terms of (R , Z) 
and ( r , z) are comparable, the compactification does not cause loss in resolution. The boundary 
conditions at r = oo or z = oo are trivial Dirichlet condition such that the metric is asymptotically 
flat. At 2 = 0, the spacetime remains symmetric about z = 0. The boundary conditions at 
r = 0 are parity conditions and local flatness condition which are consequences of axisymmetry: 
9rVtt\ r =0 = drVzn | r= o = ^rVtz\ r= Q = drVrr\ r= Q = ? ?tr| r= o = Vrz\ r= Q = IE| r _g = 9 r $| r _g = 0. 
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The simulations are performed by “standard” FDA of second order accuracy 
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(3.60) 

(3.61) 

(3.62) 

(3.63) 

(3.64) 

(3.65) 

(3.66) 

(3.67) 

(3.68) 


where the index i (or j) is to characterize the grid position in i? (or Z) direction, and At?, (or AZ) 
is the (uniform) spacing of the grids in R (or Z) direction. In the simulation, we set AR = AZ. 
The superscript n, n+l,n— 1 is to characterize the discretized time levels, while At is the spacing. 
Since compactified coordinates are used, the Courant factor is defined as Af/min(Ar) = At/AR, 
where min is taken over all A?’. The Courant factor is set to be 0.5 in the simulation. 

The residual equations are the discretized equation (3.28) (the GH formalism of GR). The 
update scheme is to obtain quantities at time level n + 1 from given quantities at level n and n— 1, 
by solving the residual equations. We solve the residual equations by pointwise Newton-Gauss- 
Seidel iteration, until residuals are smaller than a “small” threshold. 

A Kreiss-Oliger [73] style numerical dissipation is adopted to control high frequency numerical 
noises. Since we are using second order FDA, a fourth order KO dissipation is sufficient. Before 
obtaining the advanced time level n + 1, the dissipation is applied to both n and n — 1 time levels 
[61]. 
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3.4. Evolution of Massless Scalar Field under Cylindrical Coordinates 


Tests and Validation of the Numerical Scheme 

The first test we perform is to examine whether spherical initial data evolved under axisyrn- 
metric code, can remain spherically symmetric. Fig. 3.2 shows instants during the process when 
the majority of <f> field is being bounced back at the centre. The figure shows that the spherical 
symmetry is indeed preserved during the evolution. 



(a) (b) (c) (d) 


Figure 3.2: The preservation of spherical symmetry monitored by <f>. The horizontal axis is r, 
and the vertical axis is z. The simulation was carried out by 16 cores, and these graphs show 
the results at the r = z = 0 corner. Since compactified coordinates are used, the spherical 
configuration will appear to be non-spherical. However, when r and z are small, the distortion due 
to the compactification is small. This is why we choose to show the r = z = 0 corner. The four 
graphs show four instants when the pulse travels towards the center, and then gets bounced back 
and travels outwards. The graphs show that the spherical symmetry is indeed preserved during 
the evolution. 


By the general theory of numerical solutions, the numerical result is a numerical solution only 
if the independent residual tests and the convergence tests of the constraints are passed. The inde¬ 
pendent residuals are obtained by evaluating the residuals resulting from a different discretization 
which is to discretize the R and Z derivatives by forward finite difference approximation of second 
order accuracy, while the time derivative is discretized by backwards finite difference approxima¬ 
tion with second order accuracy. Fig. 3.3 shows the convergence of independent residual of the 
tt component of (3.28), and the convergence of the t component of the constraint equation. The 
results of other independent residual tests and convergence tests are similar, therefore omitted. 
These tests show that the simulation produces valid physical results, which confirms that our 
Cartesian components method does not have the regularity issue. 

Another way to justify that the solution is indeed a physical solution, is to directly evaluate 
the residuals of original Einstein’s equations R^ = kd (t^ v — without referring to the 

source functions at all. i.e. we use the generalized harmonic formalism to obtain the solutions, then 
substitute them into original Einstein’s equations to get the residuals. Fig. 3.4 shows the tt and 
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3.4. Evolution of Massless Scalar Field under Cylindrical Coordinates 




time time 



Figure 3.3: The independent residual tests and the constraint tests. Results from two resolutions 
are needed to perform the tests. Here the finer resolution’s grid spacing is one half of that in 
the coarser resolution. If the numerical results are actually numerical solutions, the independent 
residuals obtained at finer level are one-fourth of the independent residuals obtained at the coarser 
resolution, since these results are obtained from the FDA with second order accuracy, and the 
independent residuals are also discretized at second order accuracy. Fig. (a) is the independent 
residual test of the (tt) component of Einstein’s equations (3.28). Fig. (b) is the test of the t- 
component of the constraint equations H tl + T^ ~ 0. Fig. (c) is the snapshot of the (tt) component 
of the independent residual at the (r = 0, z = 0) corner, at some instant. The function with zero 
value is shown as a reference, which appears as a plane with uniform colour (light blue). The one 
taking larger value is obtained from the coarser grid, which value should be 4 times of the value 
obtained from the finer grid at all grid points at all instants. We purposely show the (r = 0, z = 0) 
corner which is the location suffering the most severe irregularity (if there is). This figure shows 
that there is no irregularity in our simulation, and the independent residual indeed convergences 
as expected at every grid point. Fig. (d) is the snapshot of the residual of constraint H t +T t at 
the (?’ = 0, ^ = 0) corner, at some instant. 


90 
























3.4. Evolution of Massless Scalar Field under Cylindrical Coordinates 


rz components respectively, as two examples. Indeed they converge as second order quantities. 

Res. of tt Comp, of Einstein Eq. 



time 

(a) 



time 

(b) 


Figure 3.4: The convergence tests using the residuals obtained from the original Einstein’s equa¬ 
tions. Fig. (a) is the convergence of the residual obtained from evaluating the (tt) component 

of the original Einstein’s equations R^ v = kg (t^ — ■ Fig. (b) is the test from the (rz) 

component. 


3.4.3 Superficially Singular Metric Representation 

To test our conjecture, instead of using (3.57), here we perform the evolution in terms of the 
following representation of metric which includes a superficially singular term 


^singular 


^ Vtt Vtr 0 

0 

0 0 rjggr 

0 0 0 

^ Vtz Vrz 0 


0 rjtz ^ 

0 Vrz 

0 0 

r/ggr 2 sin 2 6 0 

0 Vzz 


(3.69) 


where r/gg is replaced by £ • + rW . The boundary conditions of £ are f | r=0 = 0, £| r=0O = 

!> dzt \ z=0 = z\z =oc = r /( r + !)• 
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3.4. Evolution of Massless Scalar Field under Cylindrical Coordinates 


The source functions are of course 


V 




1 h t + (2/r) {rjtr/yee) ^ 

H r 


h r + (2/r 2 ) ■ £ • (r + 1 )/r]ee 

He 

= 

cot 9 

H 4, 


0 

H z j 


y h z + (2/r) {rirz/hee) / 


(3.70) 


The metric representation Singular should be problematic due to f/r terms, according to the 
conventional point of view in [50]. However, according to our conjecture, as long as £ is regular, 
and its asymptotic behaviour is compatible with the relevant FDA (in the sense that the FDA is 
able to accurately represent the derivatives of £ at small r), then the result would be regular (i.e. 
not problematic). Since T] rr ~ const + r 2 at small r, one should expect £ ~ r ■ r] rr ~ r at small r, 
and the second order finite difference approximation to its derivatives should be exact. Therefore 
the simulation in terms of £ should be regular, if our conjecture holds. 

The independent residual test and constraint convergence test of the simulation in terms of 
(3.69) and (3.70) are shown in Fig. 3.5. It shows the simulation produces valid physical results as 
well. The superficially singular term £/r does not cause problem. 




time time 

(a) (b) 


Figure 3.5: The independent residual tests and constraint tests for the simulation using the su- 
pervihally singular metric (3.69). Fig. (a) is the independent residual test of the (ft) component 
of Einstein’s equations (3.28). Fig. (b) is the test of the t-component of the constraint equations 
II^ ~ 0. These two tests show that the simulation using (3.69) and (3.70) indeed satisfies 

Einstein’s equations. 


Now we examine whether these two simulations produce the same result. We perform this test 
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3.5. Remark about Regularity 


by comparing the same quantities resulting from two simulations. Here we show the L 2 -norms of 
(a — 1) (where a is the lapse function) and <1> in Fig. 3.6, which shows that the two simulations 
indeed produce the same results. 

Therefore, our conjecture passed this challenging test. 



time 


$ from Simulations 



(a) 


(b) 


Figure 3.6: The comparison of the results from two the simulations. In the legend, “normal” means 
the result obtained from the simulation using metric representation (3.57) and “singular” means 
that obtained from (3.69) where there is a superficially singular term in the metric representation, 
(a) shows the Z/ 2 -norms of (a — 1) from the two simulations, and (b) shows the Z/ 2 -norms of <f> 
from the two simulations. These two graphs show that the two simulations indeed produces the 
same results. 


3.5 Remark about Regularity 

We have a conjecture regarding the regularity issue: provided the fundamental variables used 
in a numerical simulation are regular functions, and the FDA scheme is compatible with the 
asymptotic behaviours of the fundamental variables at the symmetric axis (or the origin), there is 
no irregularity issue. 

According to the conjecture, the key to overcome irregularity issues is to express tensors (and 
pseudo tensors) in terms of regular functions that are compatible with the FDA. One way to 
do it is to carry the regular Cartesian components into cylindrical coordinates (by coordinate 
transformation relations) and let them serve as the fundamental variables. 

There exist formalisms in the literature that are able to generate simulations without regular¬ 
ity issue, according to our conjecture. Among these specific methods, our Cartesian components 
method and Lie derivative cartoon method are equivalent (although the philosophies underlying 
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3.5. Remark about Regularity 


them are different). Both methods use Cartesian components, and are general ways to rewrite 
equations in non-Cartesian coordinates using Cartesian components, which apply to any formal¬ 
ism of GR. The background removal methods with tensorial source functions [36, 97], however, 
only apply to the formalisms where Christoffel symbols (or other connection coefficients) play 
fundamental roles, such as generalized harmonic formalism and BSSN formalism. In these for¬ 
malisms, when the background is chosen to be flat, this background removal method generates the 
same results as those of our Cartesian components method (and Lie derivative cartoon method). 
One advantage of the background removal method is that it can be applied in non-flat background, 
where our Cartesian components method (and Lie derivative cartoon method) need to be modified. 
The gauge choice in terms of the tensorial source functions needs further study. The background 
removal method used in [54, 55], however, is different from any of the above. It is also regular, 
although it is not clear how to setup coordinate gauges in this formalism. 

Above all, I emphasize the following: for the aspect of regularity, our conjecture is more 
fundamental than any specific formalisms in this chapter. For example, Cartesian components 
are not guaranteed to be singularity free—an example is the behaviour in the vicinity of the 
physical singularities. In this case the Lie derivative cartoon or Cartesian component method will 
not help. In background removal methods, the modified source function serving as fundamental 
variable, is the subtraction of two terms. In principle it is possible that both terms are singular 
but the function resulting from the subtraction is regular, and in this case the formalism can 
still generate regular results. All in all, before using any formalisms, one need to check the two 
aspects of the conjecture: regular functions being used as fundamental variables, and FDA being 
compatible with the asymptotic behaviours of the fundamental variables in the vicinity of the 
axis/origin/other special points. 
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Chapter 4 


Initial Data 


The evolution of the branewoiid is specified by the evolution equation(s) with the initial data 
and the boundary conditions. The governing equations of the braneworld are Einstein’s equations 
in the 5D bulk: R^ u — = —A g^ v . We will use the formalism of Einstein’s equations with 

the conformal factor (3.46). For the conformal factor T, we simply choose T = t/z and q = 2. 
Therefore the conformal metric is 

9a/3 4 / ^ ga/3 ( 4 - 1 ) 


According to the knowledge we obtained from the regularity discussion in Chap. 3, the most general 
yet regular metric can be 


9 = 


Vtt Vtr 
Vtr Vrr 


0 
0 

\ Vtz 


o 

o 

Veer 2 

0 

0 


0 

0 

0 

Veer 2 sin 2 ( 


Vtz 

Vrz 

0 

0 

Vzz 


(4.2) 


where fjgg is replaced by fj rr + rW. 

In this chapter and the next, we choose the unit i = fcs = 1 (which implies 87 tG 4 = 87 tG§ = 1). 
I 11 this chapter, the initial data is obtained by solving the Hamiltonian constraint, subject to the 
boundary conditions imposed by Israel’s junction condition. 

The initial data is specified on a spacelike hypersurface. Because the evolution equations are 
partial differential equations with second order time derivatives, the initial data is specified by the 
values of the evolving quantities (to be specified below) and their first order time derivatives. The 
lapse and the shift functions are related to the gauge freedom in choosing the coordinates which 
can be chosen arbitrarily, therefore the initial data is the combination of spatial (cl — l)-metric g t j , 
matter distribution, and their first order time derivatives that satisfy the Hamiltonian constraint 
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4.1. Formulation 


and momentum constraints. For time symmetric initial data, the momentum constraints are 
satisfied automatically, which leaves the Hamiltonian constraint to be solved. In most situations 
of GR, we can then choose a spatial metric ansatz with only one unknown function such as 
d l 2 = A (dr 2 + r 2 (dd 2 + sin 2 9d(j> 2 ) + d^ 2 ) . The Hamiltonian constraint is then solved for A. We 
will examine whether this situation will still be true in the braneworld. 

4.1 Formulation 

In the braneworld, Israel’s boundary condition can not be met if there is only one unknown 
variable in the spatial metric. To “absorb” this “inconsistency”, there must be at least two unknown 
variables in the spatial metric. The spatial metric can be taken to have a conformal form in (r, z) 
coordinates 

d l 2 = ^ [e 2A+2B (dr 2 + d; 2 ) + e 2A ~ 2B r 2 (dfl 2 + sin 2 0d</> 2 )] . (4.3) 

The brane is located at z = £, while z = oo is the spatial infinity. Eq. (4.3) is actually the most 
general spatial metric for the axisymmetric (in the bulk) case, as discussed in Sec. 2.4.2. The 
background spacetime (vacuum) corresponds to A = B = 0. 

The Hamiltonian constraint is 

| -1 (v’b» + ^ + 4 (i - «<«) 

+ ^_^ + | w y + § (t ,,y = o, (4.4, 

where V 2 = d rr + d zz , and U = A — B. The domain within which to solve the equations is 

(r,z) G [0, oo) x [l,oo). Let <E> be the massless scalar field that lives on brane, Israel’s junction 

condition (1.29) is then translated into 

A z = 1 - e A+B - le~ A ~ B ($ r ) 2 , (4.5a) 

D 

B, z = -^e- A - s (4>, r ) 2 , (4.5b) 

from which one see that a vanishing B is only possible when there is no matter ($ >r = 0) on the 
brane. 

The boundary conditions at the symmetric axis r = 0 are simply the parity condition and 
the local flatness condition which translate into A J n = B J n = B\ „ = < 3> r „ = 0. 

When z —> oo, the spacetime should approach the background: = B = 0. If the 
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4.1. Formulation 


matter is localized to finite r, then the spacetime should approach the background as r —* oo: 

-4|r —>-00 — B\ r — ).oo — 0 — ‘hlr-i-oo- 

There is only one equation (the Hamiltonian constraint) for two variables A and B, therefore 
there remains freedom in the initial data. We fix the freedom by imposing the form of B. Based on 
the specific forms of B, we used two specifications: Laplacian specification and direct specification. 


4.1.1 Laplacian Specification 

Since there is no other requirement on B except for the boundary conditions, we will impose an 
equation for B, for example, V 2 £> = 0. However, it turns out that there are numerical instability 
issues at spatial infinities, which suggests us to carry the asymptotic behaviour at these boundaries 
by some factors. I.e. we define Y as B = fgY where the asymptotic behaviour at spatial infinities 
is carried by /#. Since we do not know the asymptotic behaviour yet, we will try to choose /b in 
an experimental manner. There is actually a guidance we can start with. A similar problem is the 
static star configuration studied in [47], where the boundary condition at z —> 00 is asymptotically 


B ~ 


const 


(4.6) 


If the matter is localized at finite r, the boundary condition at r —> oo is asymptotically [47] 


„ const 

B ~- 

r 


(4.7) 


As discussed in Chap. 2, the asymptotic bahaviour of A at r z is supposed to be A ~ 1/rz for 
conformally flat space. We can then introduce factors like 


(r + 1 )z' 
Yr 2 

(r + z) 3 ’ 


(4.8a) 

(4.8b) 


where the r 2 in equation (4.8b) is to let the boundary conditions B | ; _ 0 = B, r | r _ 0 = 0 be satisfied 
automatically. In the numerical method, we will directly solve for X and Y instead of A and B. 
The equation of B can be imposed via Y as 


V 2 y - oY b = 0. 


(4.9) 
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4.2. Numerical Methods 


The principle part V 2 plays the role as a smoother, and the second part is to serve as an amplitude 
damper (to be explained below) when both a and b are set to be positive values, such as a = b = 3. 
The functionality of the damping term is to prevent the amplitude of Y from getting too large. 
The damping term is needed because instabilities arise at the brane when the amplitude of Y 
becomes too large. 

4.1.2 Direct Specification 

Since the only requirements on B are the boundary conditions among which only the boundary 
condition at the brane is non-trivial, we may directly set B to take the following simple form so 
that (4.5b) holds 

•/(*). (4-10) 

Z — i 

where f(z) satisfies the condition: f >z | = 1. In case we want the magnitude of B to be small, 

we can also let f(z ) satisfy a second condition: f\ z =i = 0. There is still a lot of freedom to choose 
/. e.g. we can choose it to be f(z) = (z — l)/z Pz , where p z = 4 if we want B to satisfy (4.6) (not 
necessary though). Or we can let B die off faster by choosing 


f{z) = (z- 1) • exp [ - (z - l) 2 /<t 2 2 J , (4.11) 

where a z can be chosen arbitrarily. In all the numerical results presented in this chapter, a z = 2o> 
is chosen, where ay is going to be defined in (4.14). Because (4.10) and (4.11) yield B\ z= i = 0, we 
have 

B(r,z) = ~e- A °(<h, r ) 2 -f, (4.12) 


where Ao(r, z) = A(r, z)\ z=1 . It is easy to check all the boundary conditions of B are satisfied. 

It turns out the numerical behaviour of the direct specification is better than that of the 
Laplacian specification, so we will mainly focus on the direct specification method. 


B(r, z) = 


dB 

dz 


4.2 Numerical Methods 


In order to let the spatial infinities be a part of our computational domain, compactified 
coordinates ( R , Z) are used 


R = 


r + r 0 


Z = 


-t 


z -1 + z 0 


(4.13) 
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4.3. The Numerical Solution 


where r o and zq are parameters to control the compactification. ro = zq = 1 is used for the 
numerical results that are going to be shown in this chapter. The discretized grids are uniform 
in R and Z , and the equations are discretized by finite difference approximation with second 
order central stencil, which are shown by eq. (3.60) to eq. (3.68). Then Newton-Gauss method is 
employed iteratively to solve for the numerical solution. After the numerical solutions are obtained, 
we validate the solutions by the mean of independent residual test. 


4.3 The Numerical Solution 


4.3.1 The Solution and Apparent Horizon 

We used the following profile for the matter field 


<f> = stf ■ exp 


(r - Zq) 2 ] 


(4.14) 


The condition <I>„ = 0 needs to be satisfied in the initial data, therefore we either choose 

>' lr=0 

xo = 0, or choose xo to be at least several ay. 

The independent residual, calculated as the residual of equation (4.4) under a different dis¬ 
cretization scheme (other than eq. (3.60) to eq. (3.68)) that is of second order accuracy, should 
behave like a second order quantity everywhere in the domain if a numerical solution is obtained. 
To show the independent residuals compactly, we use their L 2 norms. As an example, we per¬ 
form the independent residual test to the numerical results obtained from the calculation with 
(&/, cr r , xo) = (2.8,0.15,0) using the direct specification. The L 2 norms of the independent resid¬ 
ual at resolutions 64 x 64 (the domain in R direction and Z direction are uniformly divided into 
64 intervals), 128 x 128 and 256 x 256, are 0.0799,0.0194 and 0.00498, respectively. The L 2 norm 
shrinks by a factor of 3.9-4.1 when the grid spacing decreases by a factor of 2. The independent 
residual indeed converges to zero at second order, therefore the numerical scheme is justified. 

The results of the numerical calculation with {srf,<j r ,x 0 ) = (3.2,0.3,0) using the direct speci¬ 
fication, are shown in Fig. 4.1. The scalar field is strong enough to produce an apparent horizon 
which is shown in Fig. 4.2. 


4.3.2 Brane Geometry as Seen by a Brane Observer 

In this subsection we focus on the geometry of the brane, since the bulk is not directly observ¬ 
able. In the next subsection, we will study the geometry of the bulk. Again, here we reiterate 
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4.3. The Numerical Solution 


Function A Function —B 



r/(l + r) r/(l + r) 

(a) function A (b) function —B 


Figure 4.1: Function A and — B in the initial data metric (4.3) obtained from the parameters 
[pf > cr r , xq) = (3.2,0.3,0) 


Apparent Horizon 



Figure 4.2: The apparent horizons for the configuration with parameters (&/, cr r , Xo) = (3.2, 0.3, 0). 
The red x in the figure is the apparent horizon calculated using the brane geometry only. It does 
not coincide with the intersect of the brane and the bulk apparent horizon. 
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4.3. The Numerical Solution 


that the masses of the spacetime obtained from the brane geometry, play no direct roles in the 
braneworld. Especially, they do not really represent the masses of the brane. The purpose of 
studying these masses, is to compare with GR and obtain some observable difference from GR. 

We will study two different masses defined on the brane—brane ADM mass and the Hawking 
mass (see Sec. 2.5.1 and 2.5.2). The proper area radius r on the brane is 

f = e A ~ B r, (4-15) 


so that the metric on the brane is now 

dSbrane = 9ffAf 2 + f 2 dff 2 , (4.16) 

where g?? = e 4B [1 + r (AA/Ar — AB / Ar)\~ 2 . 

In terms of gff and f, the Hawking mass is 

Mu • ^T-(l (Srr) -1 ) , (4.17) 

which has f dependence. Then the r dependence of the Hawking mass is an observable difference 
from 4D GR. One example is shown in Fig. 4.3. 

The Hawking mass goes to ADM mass as f —> oo 

Madm = _lim (l - (ffff) -i ) = ni -t " jl ■ (4.18) 

and /?i are defined via the asymptotic behaviour of A and B at large r, which are assumed to 

be A |_« ot\ / r and B | _ « /?i/r. As it is shown in Sec. 4.3.3, the asymptotic behaviour of B 

actually implies f3\ = 0. 

The masses for the results of the configuration with (&/, cr r , xo) = (3.2, 0.3,0), are shown in 
Fig. 4.3. In the next section we will study the total mass in the bulk. To distinguish different 
masses, the ADM mass on brane will be referred as Mb ra neADM- 

4.3.3 Asymptotic Behaviour 

The asymptotic behaviour at r z is crucial for the calculation of total energy (see Sec. 2.4.2). 
The construction of B in the direct specification according to eq. (4.12), directly implies |B| |A| 

at the r z region. When B is negligible, in Sec. 2.4.2 we showed that a solution at r z region 
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4.3. The Numerical Solution 


Masses on Brane 



■■-■Mb 

' Midm 

. Madm 

— apparent horizon 


0.4 


0.6 


0.8 


l+r 


Figure 4.3: Brane Masses from the initial data configuration with parameters (&/, a r , Xo) = 
(3.2,0.3,0). The result is Madm = 22.67. The graph shows that the Hawking mass is not a 
constant in the braneworld, in contrast to 4D GR. Hawking mass agrees with ADM mass only at 
r —> oo. Madm (defined in equation (2.58)) still blows up around the apparent horizon. 
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4.3. The Numerical Solution 


was 


A ss — when r z. (4-19) 

rz 

However, as discussed in that section, there is no proof of the uniqueness of this solution. Therefore 
in this section we need to test whether the solution is indeed given by eq. (4.19). Please refer to 
Fig. 4.4 for the results of a series of simulations from the family (ji/,cr r ,x o) = (. 2 /, 0.3,0) using 
the direct specification. The graphs show that the asymptotic behaviour is indeed described by 
eq. (4.19). 



10° io 1 io 2 



10° io 1 io 2 






Figure 4.4: The asymptotic behaviour of A at r z region, shown by the simulations from the 
family (&/, <j r ,xo) = (.e/, 0.3,0). In the first diagram, we plot r x A versus z in log-log scale. The 
red line, is plotted assuming r x A = et\/z. If A m a\/rz , one would expect the plots for different 
r should follow the red line. Other diagrams are plotted in the same way, but with different a\ 
resulted from different stf. The diagrams look almost identical, despite that a\ has changed over 
three orders of magnitude. It means that the asymptotic behaviour of A at r z region is indeed 
A ~ ai/(rz). 
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4.4. The Total Energy 


4.4 The Total Energy 

The asymptotic behaviour of eq. (4.19) is confirmed by Fig. 4.4, therefore the calculation in 
Sec. 2.4.2 applies. The calculation shows that the total energy (M to tai) in the whole spacetime is 

~ 8-o,. (4.20) 

Consequently, as shown by eq. (2.67), the relation between the brane “energy” and the total energy 
is then simply 

A^braneADM — A/total- (4.21) 


4.4.1 The Relation with the Area of Apparent Horizon 

In 4D GR case, there is a simple relation between the Hawking energy of the spherically 
symmetric vacuum spacetime and the area of the horizon of the black hole sitting at the symmetric 
center, which is eq. (2.62). In the braneworld, we would like to examine whether the total energy 
can also be characterized by the area of apparent horizon in 5D braneworld. To study this relation, 
we choose the configurations where the matter outside of the apparent horizon are negligible. 

Fig. 4.5 shows the relation between the total energy M to tai> and \/47rAbuik, where Abuik is the 
area of the apparent horizon in the bulk. The figure shows these two quantities are equal. Therefore 
we have 

Altotai = v/TrAbuik- (4.22) 

How to understand this relation? Eq. (4.21) tells us that the mass can also be understood as 
the ADM mass calculated on the brane. The size of BH—the areal radius of the intersect of the 
horizon with the brane—is r a ~ (0.7, 2.0) according to the area-versus-radius relation shown by 
Fig. 1.1, where r a is the areal radius of the horizon on the brane. These BHs are then considered 
to be medium to large. For large BHs, the area of the horizon in the bulk reduces to the area 
of black string, which is equal to the area of the horizon on the brane, and the brane geometry 
(without matter) reduces to 4D GR. Eq. (4.22) might be just 

AfbraneADM = \jTu Ab rane* (4.23) 

On the other hand, given there are many BHs of medium size in this data set, the knowledge 
of large BHs (that the areas reduces to those for black strings) might not apply. Therefore the 
relation might be generically about total energy and the area of apparent horizon in the bulk. 
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4.4. The Total Energy 


Mass-Area Relation 



Figure 4.5: The mass-area relation of apparent horizon. This plot is produced from simulations 
with various parameters (£/,(T r ,xo), using both direct method and Laplacian method. The only 
criteria is that the configurations defined by the parameters are able to produce apparent horizons 
in the bulk. 
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4.5. Discussion: the Results of Laplacian Specification 


At the end, we emphasize that the range of the size of the BHs is merely r a ~ (0.7, 2.0), and 
Fig. 4.5 exhibits a small but systematic deviation between the data and the line of a linear fit. 
Therefore further study is needed. 


4.5 Discussion: the Results of Laplacian Specification 


The total energy is obtained in eq. (4.20), which was derived based on the condition that 
|B| -C \A\ in r 2 region. By the construction of the “direct specifiction”, \B\ -C |A| is met 
automatically. On the other hand, without pre-setting this condition, this condition emerged in 
the star solution [47] and the small black hole solution [14]. In our initial value problem, there is 
no such a pre-set condition in Laplacian specification method. It is then interesting to see whether 
this condition can still emerge. Fig. 4.6 shows the results of a series of simulations from the family 
(.s/, ay, a?o) = {sf, 0.15,0.5) using Laplacian specification, where we see that indeed A « a\/rz 
emerges. The validity of the discussion leading to A/total = ay/G 5 is based on the assumption 
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Figure 4.6: The asymptotic behaviour of 4 at r > 2 region obtained via Laplacian Specihcation, 
shown by simulations from the family = ( 22 /, 0.15, 0.5). This diagram shows that the 

asymptotic behaviour of A is indeed A ss a\/{rz). Instead of almost identical configurations as 
in Fig. 4.4, this diagram shows interesting pattern: 1) the asymptotic behaviour of stronger data 
is more like A ss a\/rz\ 2) disregard the trend with the strength of the data, as long as r is big 
enough, there is always A ss a\/rz. The asymptotic behaviour of B is observed to be B ss c/r 23 . 
The fact that c/a 1 decreases with ai explains that A become more like a\/(rz) with the increase 
of ay. 


A « ay /rz and |£?| -C |A|. Since B ss c/r 2 3 and A ss ay /rz, V (small) e > 0, 3 Q = [cv/a ie ) 10//3 
such that | B/A | < e for r = Q,z G [1 ,v ■ Q\ (which is the region to calculate the total energy, 
as shown in Fig. 2.2(a) and discussed in Sec. 2.4.2). Therefore the assumptions hold, and we still 
have Mtotai = ay/G 5 - 
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4.6 Prepare the Initial Data for Evolution 

Often the lapse a and the shifts fju are freely specifiable in initial data. However, for the 
braneworld, the Israel’s junction condition imposes non-trivial boundary conditions on the lapse 
function and the shift functions. For our specific initial data ansatz with time symmetric initial 
configuration, the boundary conditions for the shift functions are automatically satisfied by setting 
fjti = 0. The lapse function, on the other hand, has a non-trivial boundary condition. We denote 
the tt component of the metric as gu = — (L 2 /z 2 ) e 2a , then Israel’s boundary condition for a is 

a z = 1 — e A+B + ^e~ A ~ B ($, r ) 2 . (4.24) 

Using the direct specification, a can take 

« = /•( 1 — e A + ^e - 4 ($ , r ) 2 ) , (4.25) 

where / is the same / used in eq. (4.10) to specify B. 

4.6.1 Prepare the Initial Data for Specific Initial Source Functions 

Given the values of the metric functions and their time derivatives, the constraints for source 
functions h p ~ 0 (which is satisfied at initial instant), can give us the initial values of the 

source functions. Alternatively, as used by Pretorius [63], if one would like to set a certain gauge 
at the initial instant by setting the value of then one can use the constraint equations to work 
backwards to get the values of dtfjty.. In this way the initial values of dtfjtfi is expressed in terms of 
the value of the metric components and their derivatives. Let us denote the obtained expression 
as fjt^tiv, dfj). 

This procedure works well in ordinary spacetime (such as the 5D massless scalar field studied 
in section 3.4). However, in the braneworld, generally dfj) is not continuous at the brane 

under the perpendicular gauge (2.7). This is because, dtfjtz should be zero at the brane, due to 

the boundary condition fjt z \ _ 1 = 0 at all time. On the other hand, fjtz,t(fj,dfj) will include z 

derivatives of other metric functions (i.e. fjtt,z,Vrr,z, fizz,z, tide,z), which generally results in a non¬ 
zero fjtz,t(,V,dfj) at the brane, due to Israel’s junction condition, i.e. fjt z ,t is not continuous at the 
brane. 

This discontinuity problem can be solved by properly adjusting fj zz . If we do not choose 
fizz — fjrr at the initial instant, there is the freedom to setup the value of ij zz since neither Israel’s 
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junction condition nor parity condition impose any requirement on fj zz . We can use this freedom 
to adjust fj zz such that f/tz.tiv, dfj) = 0 at the brane. In this case the simplest ansatz for the initial 
metric is 

ds 2 = ^ — e 2tt df 2 + e 2A+2B dr 2 + e 2A ~ 2B (d(9 2 + sin 2 9d<j> 2 ) + e 2b dz 2 ^j , (4.26) 


for which the Hamiltonian constraint reads 
_l +e 4B 3e 2(A+B) 1 


2 r 2 


+ J^ e2{A ~ b+B) ( 6 + 6 * 2 (A.) 2 + 3 zB tZ + 2z 2 (B z ) 2 


“h zb^ z (3 4- zB^ z ) — zA, z (9 H- 3zb z “l - 4 zB z) ~\~ 3 z^ A ^ Z z — ^B z 

1 . . 2 1 1 / . 2 2^4 tr -|- b r — 4 B r 

+ 2 (A**) 4- 2^’ r ^’ r 2 - 1 ---~ 

3 5 1 

— 3A r B r — -6 r £>. r — (B t r) + A, rr -f- ~b,rr — B rr — 0. 


(4.27) 


The boundary conditions of a, A, B are now 






-2A-2B 


($,r 


B - = -r. 


-2A-2B 


($,r 


(4.28) 


The only requirement on b (from solving fjt Zt t(fj, dfj) = 0) is its derivative at the brane 

b. z = - (h z + l e ~ 2A ~ 2B ($, 2 ) 2 ) , (4.29) 

so that dtfjtz is continuous on the brane. The h z in (4.29) is the target initial value of h z that we 
try to achieve at the initial instant. For example, h z = 0 if the harmonic gauge is chosen to be 
imposed at the initial instant. 

Similar to the direct specification method applied to B above, we can adopt direct specification 
method to specify (a, 5, b) 

« = /• (^ e “ 2A (*T) , B = f ■ (-\e~ 2A (4>, r ) 2 ) -W-( (/Vi ie“ 2A (<b,,) 2 )) . (4.30) 
Or 

a = f ■ (<M 2 ) ,Bmf- (Te - 2 - 40 (<f>, r ) 2 ) , b = f • (- (h z + (* ;) 2 )) , 

(4.31) 

where Aq = A\ z= i, is the value of A on the brane. While both of (4.30) and (4.31) work, the 
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numerical experiments show that (4.31) enables the Hamiltonian constraint equation (which is an 
elliptic equation in term of A) converges faster when performing numerical relaxation. Therefore 
we adopt (4.31). 

For an example of the initial data, please refer to Fig. 4.7. 


Function A 



0 0.2 0.4 0.6 0.8 1 

r/(1 + r) 


Function —B 



r/(l + r ) 


(a) Initial Data: function A 


(b) Initial Data: function —B 


Figure 4.7: Function A and —B in initial data metric (4.26) from the initial data configuration 
with parameters (<£/, <r r , xq) — (2.5,0 .2,0) 


4.6.2 Total Energy 

The total energy from the metric (4.26) can be obtained by the same method presented in 
Sec. 2.4.2. For this metric, according to equation (4.30) or (4.31), we have b ~ 0 and B ~ 0 at 
large r (or large z) region. Therefore the embedding conditions (2.40) and (2.41) reduce to 


r 


z 

( f ') 2 + ( z ') 2 

~2 



(4.32a) 

(4.32b) 


where ' = d/dz. 

Only the asymptotic behaviour of A at r A> z region contributes to the calculation of total 
energy. The asymptotic behaviour is (shown by Fig. 4.8) 

A « — at r>z. (4.33) 

r 
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This asymptotic behaviour gives the solution to (4.32) at the lowest order of A as 


z 


z, r 


r • exp 



which yield the total energy 


Mtotai = CX 1 /G 5 = 87rai. 


(4.34) 


(4.35) 


Therefore, the equality of the total energy with the ADM “mass” calculated on the brane, also 
holds for metric (4.26). 




Figure 4.8: The asymptotic behaviour of 4 at r > z region, shown by numerical solutions for 
the family (^/,cr r , xq) = (&/, 0.15, 0.5). The diagrams are almost identical, and here we only show 
the one with the smallest aq and the one with the largest aq. These diagrams show that the 
asymptotic behaviour of A is A « aq/r for a\ over a wide range of magnitude. Therefore, the 
conclusion of A « 21 at r z, is robust. 
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Chapter 5 


Evolution 

5.1 The Evolution as an Initial Value Problem 


Axis 



AdS 


Brane [z = l\. Israel’s BC and matter’s equations of motion. 


Figure 5.1: The evolution as an initial value problem with boundary conditions, z G [l, oc) is 
the extra dimension, while the brane is located at z = t (expressed by the shadowed line in the 
diagram), r G [0, oo) is the radius coordinate, and r = 0 (z axis) is the symmetry axis of the SO(3) 
symmetry. 


The evolution of the braneworld is specified by the evolution equation(s) with initial data and 
the boundary conditions. The initial data was discussed in Chap. 4. Here we present the evolution 
equation and the boundary conditions. 

The governing equations are composed of the 5D Einstein’s equations in the bulk and the 
matter’s equation of motion on the brane. The 5D Einstein’s equations are 


Gaft — ^9cx/3 • kdJ'oifi 


Rql(3 — Qoc/3 


T 


d- 2 


(5.1) 
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Its form in terms of the conformal metric g a p is 


R(x(3 — fod, [ ^ot.fi Qot.fi 


T 


d- 2 


- R 


(») 

otfi ’ 


(5.2) 


where g , 4/, R^p and R a p were introduced in section 3.3. For the conformal function and the 
conformal parameter, we take q = 2 and T = 1/z in our calculations. 

The most general conformal metric in terms of regular functions ?y, is (see eq. (3.35)) 


( 


Vtt 


Vtr 


9a/3 — 


Vtr Tjrr 

0 0 

0 0 


y fjtz Vrz 


o 

o 

Veer 2 

0 

0 


0 fjtz ^ 

0 Vrz 

0 0 

Veer 2 sin 2 6 0 

0 fj zz I 


(5.3) 


where fjgg is taken as equal to f) rr +rW , which makes the local flatness condition vee\ r _ Q = Vrr | r _ 0 
be satisfied automatically. 

The evolution equation (5.2) can be solved in two ways. The first way, is to use (3.48), in terms 
of the following source functions (see eq. (3.36)) 




^ ht + (2/r) {fjtr/veg) ^ 

H r 


h x + (2/r) (Vrr/vee) 

He 

= 

cot 9 

h 4> 


0 



y h z + (2/r) (vrz/vee) ) 


Alternatively, a second way, is to use the generalized harmonic formalism with conformal function 
and tensorial source functions, which is equation (3.52) 


- ^V 0 V,)V - V Q ^( Al V^ ) ff Q/3 - V ( Ji v) + h a C a ^ - C a ^& au 
• k d (t^ - - U,,,, U { *\ (5.5) 

When the background in (5.5) is taken to be flat spacetime (which is what we adopted in the 
numerical calculations presented in this chapter), the two ways are equivalent. For discussion and 
presentation purpose, we will adopt the second, eq. (5.5). 
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A second part of the governing equations is the equation of motion of the matter on the brane, 
which only “feels” the 4D metric. For massless scalar field, the equation of motion is 


= 0 , 


(5.6) 


where D is the covariant derivative associated with the brane metric h a p. 

Next we look at boundary conditions. The boundary conditions at the spatial infinities (r —* oo 
or z —> oo) are Dirichlet conditions such that the metric takes the form of the background 

ds 2 = 3 ( - d t 2 + dr 2 + r 2 (dd 2 + sin 2 Odcf) + d;r 2 ), (5.7) 


which is a Poincare patch of the AdS spacetime (see Sec. 1.4). The background of RSII braneworld 
is the z > t portion. In terms of the conformal metric components, the boundary conditions at 
r —» oo and z — » oo are fj rr = 1, fjtt = — 1, rj zz = 1, fjtr = 0, rjt z = 0, fj rz = 0, W = 0, and <f> = 0 for 
the matter at r —> oo. 

The boundary conditions on the symmetry axis are parity conditions, and the local flatness 
condition, which are | r=0 = 0 , fj rr>r \ r=0 = 0 , fj tz ,r\ r=0 = 0 , ?% z , r | r=0 = 0 , fj tr \ r=0 = 0 , fj rz | r=0 = 
O.U'ir 0 = 0, T,r| r=0 = o. 

The boundary conditions at the brane are Israel’s junction condition (1.29) 


^ = -k d (A 


a/3 


^a/3 

'd ~2 


4“ Fa/3 h a p 


d -2 


(5.8) 


which translates into conditions on ija. z , ijtr.z , Vrr.z ? . The expressions are long and the specific 

forms do not matter at this stage, and are thus not written here. The expressions are in terms of 
metric functions, their first order derivatives with respect to r and t, and 4> ir , <3?^. As discussed in 
Sec. 2 . 1 , generically there is no boundary condition for fj IIZ , yet we need their boundary conditions 
in the numerical calculation. These conditions are more subtle and will be discussed below. 

The massless scalar field <f> only exists at z = £. Therefore it is not defined in the bulk. 
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5.2 Details of the Numerical Implementation 

For spatial coordinates, we use “compactified coordinates” to include spatial infinities into the 
calculation domain, and numerical grid is uniform in R and Z: 


R = 


r 

r + r 0 ’ 


Z = 


z-i 

Z - I + Zq’ 


(5.9) 


where r o and zq are parameters to control the scale of the compactihcation. The derivatives in 
the equations need to be changed accordingly. For example, /,r = f fo • § 7 , etc. The numerical 
calculations are performed using the central stencils of finite difference approximation operators 
of second order accuracy, whose explicit forms are 


dtf -t '*’ J 

dnf 

dzf -» 

dttf —> 
^Rflf 


fn -\-1 rn 

J i j J i ■ 


®zzf 
d m f 

d tzf 

dfizf 


2- At 

fn _ fn 

Ji+1,3 

2 • A R 

fn _ fn 

J i,j +1 

2-A Z ’ 

i*n+l _ 9 fn 1 fn— 1 

J i,j "«/ i,j ' J i,j 

lW 2 ’ 

fn _ 9 fn fn 

Ji+l,j ' Ji—l,j 

H 2 

fn _ 9 fn fn 

J' Ji,j — 1 

(AZ 


4 • At • Ah’ 

1 

4 • At • AZ 
1 

4 • Ai? • AZ 


/ /»n+1 tn+1 _i_ fn— 1 fn— 1 \ 

Vyi+ljj ‘ Ji— l,j 


( fn -\-1 /*n+l , fn— 1 fn— 1 \ 

/rn fn , fn fn \ 

\Ji+l,j+l Ji-lj+1 “T /i-lj-l Ji+l,j—l) • 


(5.10) 

(5.11) 

(5.12) 

(5.13) 

(5.14) 

(5.15) 

(5.16) 

(5.17) 

(5.18) 


where the index i and j are grid indices in R and Z directions, and A R and AZ are the (uniform) 
spacing of the grids in R and Z directions. The superscripts n, n + 1, n — 1 are the discretized time 
levels, where At is the spacing. In the simulations performed in this chapter, we set ?’o = z$ and 
the resolution in R is the same as that in Z. Courant factor Af/min(Ar) = Af/(ro • AR) is set to 
be 0.5. 

The residual equations are the discretized eq. (5.5) and eq. (5.6). Let n denote the current time 
level. The update scheme is to obtain time level n +1 from given quantities at level n and n— 1, by 
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solving the residual equations. We solve the residual equations by pointwise Newton-Gauss-Seidel 
iteration in a black-red manner (see, e.g. [73]) until residuals are smaller than a “small” threshold. 

A Kreiss-Oliger (KO) [73] style numerical dissipation is added to control high frequency numer¬ 
ical noises. Since we use second order FDA to discretize eq. (5.5) and eq. (5.6), a fourth order KO 
dissipation (see [61] for the specific form) is employed. Following [61], the dissipation is applied to 
both n and n — 1 time levels before solving for the advanced time level n + 1. 

Adaptive Mesh Refinement (AMR) is used to reach high resolution (when needed). We used 
PAMR/AMRD [65] as the tool to realize the parallelization of the code. Both of the KO style 
dissipation and AMR are built into PAMR/AMRD. The simulations producing BHs in this chapter, 
however, were obtained using only one CPU with uniform grid structure, because a shooting 
method is employed to locate apparent horizons, which is not parallelized, nor adapted to AMR. 
Our plan for the future is to upgrade the code to use flow method to locate apparent horizons, 
which can be parellelized and adapted to AMR [62]. 

5.3 Gauge Freedom 

As discussed in section 3.3, there are two ways to perform evolutions using GH formalism: the 
BSSN-like method and the source function driven method. Here we adopt the source function 
driven method. Therefore we need to consider how to impose gauges via the source functions, 
and how to make sure the + AP M ~ 0 constraints are preserved during the evolutions. In this 
section we focus on imposing gauges. Please refer to section 5.4 for the study of the constraint 
preservation. 

5.3.1 Fixing Gauges via Source Functions 

There are gauge freedoms in gravitational theories, which are about the coordinate choices. It 
is important to properly specify coordinates in numerical relativity to avoid coordinate pathology, 
to evolve spacetime with strong fields, and to deal with physical singularities. For the generalized 
harmonic formalism, the following gauges are generally adopted in the simulations in the literature. 

(1) The simplest gauge is the harmonic gauge 


= 0 . 


(5.19) 
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(2) The lapse driving gauge used by Pretorius [63]: 

n~h t + d - c 2 ht^ = 0, (5.20) 

a s 

where s > 0, Ci > 0, c 2 > 0. It is generalized from the damped harmonic equation x^t + Ci • x + 
c 2 • x,t = 0. With □, the equation has the functionality as a smoother. The effect of this gauge 
is to drive the lapse function towards its desired value do- For example, Pretorius found that 
the instability tends to happen at the apparent horizon excision boundary when the value of 
the lapse function is too small. He then chose do = 1, the value of the lapse function in flat 
spacetime. 

One can also apply this method to the spatial components hi, to achieve the desired gauges. 

(3) Damped wave gauge [60, 98]. In our case, the gauge reads 

h M = ci log (^-^j - c 2 Vuifr/a, (5-21) 

where fj = — (?) rz ) 2 J ( fjgg) 2 is the determinant of the spatial (conformal) metric in 

Cartesian coordinate. P, c\, c 2 are positive parameters. The effect is to damp out the dynamics 
in spatial coordinates on the time scale l/c 2 , and to suppress the growth in \fr\ld (when 
P = 1/2) on time scale 1/ci- 

In ours simulations, (3) was adopted for the simulations that produce small black holes. For 
small black holes, it is crucial to set both r$ and zo much less than 1 to let the black hole boundaries 
include many grid points, so that the small black holes are well-resolved. (2) was adopted for 
the simulations that produce large black holes. Pretorius’ group successfully performed many 
simulations [62, 66, 67] using (5.20) as the slicing condition, and hi = 0 for spatial coordinates, 
which are also the conditions we use. do = 1 is chosen. For large BHs, (1) also works well. In any 
case, after BHs are obtained and the evolution is stablized, we gradually change the gauge into (3) 
to drive the value of d towards do = 1 so that the lapse rate of the coordinate time is comparable 
to that of the proper time, which enables us to define and to monitor the apparently stationary 
state that is going to be introduced in Sec. 5.7.1. 
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5.4 Constraint Violation and the Cure 


5.4.1 Constraint Damping 

The generalized harmonic formalism, has the constraints = H^+T^ — 0 (or C = h^ + AT^ — 
0). During the evolutions, there are always numerically errors that violate the constraints, and 
often the modes of the deviation from the constraint equations grow with time, and drive the 
evolution away from Einstein’s equations. The phenomena, numerical violation from constraint 
equations growing with time, is a very common challenge for numerical relativity (see, e.g. [58]). 

One way to improve the performance, is to add the following constraint damping terms [58] to 
the left hand side of the evolution equations (3.28) 

Znv = k (n(u,C v ) - ’—^-gnv' n PCfj \ , (5.22) 

where k > 0 and p > — 1 are constant parameters. For evolution equations (5.5), the damping 
terms are changed to 

Znu = k (n^C v) - ^-^-9nun p C^\ , (5.23) 

where + AT M . These damping terms have been proven to be useful in keeping the 

constraints satisfied for simulations in “ordinary” spacetimes (e.g. the simulation in 5D carried out 
in Sec. 3.4). However, for the braneworld, it turns out that these damping terms are not sufficient 
— the constraints do not converge to zero as the resolution increases. In fact, the residuals of the 
constraints almost stay the same when the resolution increases. 

We experimented with many changes to the damping term—some worked better than others. 
One version was to make k spacetime dependent. We tried a few choices and it turns out the 
following choice worked to a certain degree 


K —K 


(z — C () r 


(5.24) 


where c was chosen to be, for example, 0.9 or 0.95, which effectively puts more damping near 
the brane. This change enhanced (a little) the convergence of the constraints. However, the 
convergence criterion is still not perfectly satisfied. 

It turns out the following choice is more successful, although we can not explain why at this 
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moment. More work is needed. 


Zi_lv — ft 


n^C u ) - F ■ 


(5.25) 


where F may depend on spacetime, i.e. only the second term of Z M „ is multiplied by F , which 
is a spacetime dependent function. We experimented with different forms of F , and the following 
family worked the best so far (among all our trials) 


F = 


z n — c £ n 


(5.26) 


where c is close to 1, which again adds more damping on the brane. A choice was n = 4, c = 0.99 
(refer to Fig. 5.2). This damping was successful, since the constraints converged at the expected 
order. 

The success of this method means that the second term in constraint (5.23) and the first term 
may have very different effects. 

The constraint damping is very important. We do not have specific guidance to give the forms 
of such damping parameter settings. A survey over the parameters showed that the effect of 
damping is very sensitive to parameters. Also, even for the most successful result (Fig. 5.2), the 
independent residuals (and the residual of constraint equations) increase rapidly after certain time, 
i.e. very long term simulations seems to be problematic. 


5.4.2 Imposing Constraints on the Brane 

From the result of Fig. 5.2, we learn that 1) the gauge choice /t ; , = 0 works; 2) the damping 
should be larger at the brane, which motivates us to exactly impose constraints /t /t + = 0 on 

the brane. 

How to impose constraints h fl + AF M = 0? What one usually does in successful simulations in 
ordinary spacetimes (non-braneworld, such as the 5D simulation in Sec. 3.4), is to set h's to certain 
pre-set values (e.g. in the case of the harmonic gauges, we have h = 0). In this situation, the effect 
of constraint damping is to drive the metric components to the values that satisfy the constraint 
equations, rather than doing something to the source functions h^. i.e. to impose constraints, it 
is the metric (rather than h's) that should be guided. 

How does this guidance happen? Let us define should /i M = — a p — F" a/3 ^ 9 0,3 ■ he. they 

are what the h's should be, if the constraints are exactly satisfied. The should /ps are functions of 
the metric components and their derivatives. One can set the value of metric components or their 
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Res. of tt Comp, of Einstein Eq. 




(a) (b) 





(e) 


Figure 5.2: This convergence test is for the results obtained from the simulation using harmonic 
gauge = 0, and constraint damping Z^ v = n{n^C v ^ — ^ E g IJ vh a C a ■ —rrogg)• The initial data 

is $ = exp [— (r — xq) 2 /<j^ with (sV, ro,oy) = (0.05,1,0.25). After the results are obtained 
by the generalized harmonic formalism with a conformal function, the results are substituted into 
the original Einstein’s equations in terms of the original metric functions without the conformal 
function, to get the residuals. The residuals should converge to zero as second order quantities, 
if the results are numerical solutions. Fig. (a) is the convergence of the residual obtained from 

the (tt) component of Einstein’s equations R = kd (t^ — vjrjj. Fig. (b) is the test from 

the (rz) component. Fig. (c,d,e) show the convergence of constraints: h M + AF^, with p = t,r,z, 
respectively. The residuals converge at the expected order. However, this is the result obtained 
from very extreme damping, and all the residuals have up-climbing tails. 
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derivatives at the brane to let should /i M = h M (therefore imposing constraints), which can be done 
in multiple ways. Here we adopt the following. shoulc ri M = can be equivalently expressed by 
setting the values of fjtz,z,Vrz,z and fjzz,z, which in principle can serve as the boundary conditions 
for fj liz —as shown in Sec. 2.1, generically there is no boundary condition for rj^ z . 

In addition, we impose the perpendicular gauge at the brane ( fjt z |_ e = fjrz |_ e = 0) since this 

gauge gives the smoothness of the apparent horizon across the brane, as discussed in Sec. 2.2.2. 
Now we have two conditions for fjtz (and two conditions for fj rz ): one is constraint imposing 
condition which is a condition on the value of fjt z , z (and fj rz ,z)i the other is the perpendicular 
gauge condition which is a condition on the value of f] tz (and rjrz)- Imposing two boundary 
conditions on one function is achieved by the following trick: let one of the conditions be satisfied 
automatically by choosing the forms of the functions: 


fjtz = (z - i) ■ it, (5.27) 

fjrz = (z-t)- ir- (5.28) 

In this way, the ijtz\ z= e = fjrz \ z= c = 0 are automatically satisfied. The conditions on fjt z , z and 
fjrz.z are now converted into the conditions on the values of and £ r . We then use £ t and as 
fundamental variables instead of rj± z and fj rz . 

It turns out that this method works well, and there is no need to use special damping. 
Beyond the constraints imposed on the brane, the “normal and plain” damping term v = 
k [n^C v ) — ^r9iivfi a Ca'j is still employed to control “ordinary” violation modes in the bulk. The 
tests of the method are shown in Fig. 5.3. 

5.5 The Evolution with an Apparent Horizon 

During the evolution, sometimes singularities are formed and the code crashes. There are at 
least two types of singularities. The first is the coordinate singularity due to pathological coordinate 
gauges, which can be avoided by properly choosing coordinate gauges (which is non-trivial). The 
second is the physical singularity. However, Penrose’s cosmic censorship hypothesis [84] states 
that, the physical singularity is always hidden behind an event horizon. While there is no proof of 
this hypothesis, it does seem to be satisfied in many cases. 

Event horizon is the boundary in spacetimes that separates the events which can be causally 
connected with future Poincare horizon by future oriented null geodesics, from the events which can 
not. The interior surrounded by event horizons, by definition, can not affect the regions outside of 
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(a) 


(b) 



time 


Convergence of h T + Ar r 



time 



time 


(e) 


Figure 5.3: The convergence Tests for the simulation using gauge = 0. The initial data is 
$ = exp [— (r — a;o) 2 /cr^] with (j^/,ro,a r ) = (0.05,1,0.25). After the results are obtained by 
the generalized harmonic formalism with a conformal function, the results are substituted into 
the original Einstein’s equations in terms of the original metric functions without the conformal 
function, to get the residuals. The residuals should converge to zero as second order quantities, 
if the results are numerical solutions. Fig. (a) is the convergence of the residual from the (tt) 

component of Einstein’s equations R^ v = kd ~ 9^ Fig. (b) is the test for the (rz) 

component. Fig. (c,d,e) show the convergence of constraints: h^+AT^ with p = t, r, z, respectively. 
In all the tests shown in the figures, the spacing of the coarser grid is A R = A Z = 1/256 and the 
spacing of the finer grid is A R = A Z = 1/512. The simulations were performed using 16 CPUs 
and the test is performed on the result obtained by the CPU at the R = Z = 0 corner, which is 
the region that would suffer from the most severe problems (if there was). The convergences are 
shown to be good. After t ~ 2, the residuals suddenly decreased, because the interesting dynamics 
propagated out of the region where we evaluate the residuals. 


121 
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the event horizons. Therefore, one way to avoid the physical singularity in the calculation domain 
is to perform the evolution without referring to the interior of the event horizons (so that the 
interior of the event horizons is excised from the calculation domain). 

Event horizon, however, is a concept based on the global structure of the whole spacetime, 
which is therefore not quite useful during the evolution since one can not tell the event horizon 
until the full evolution is completed. But the full evolution is not able to be obtained without 
knowing the event horizon to excise the physical singularities. Fortunately there is the concept of 
apparent horizon which is locally (in time) defined. Apparent horizons are not guaranteed to exist 
in a certain evolution. But, if they do exist, they are inside of the event horizons. Therefore the 
interior of an apparent horizon can not affect the exterior of the event horizon. Since the apparent 
horizon lies inside the event horizon, and sometimes by a long way, the Penrose hypothesis does 
not guarantee that singularities can not exist outside the apparent horizon. But again, often 
the physical singularities (to be formed in a future instant) are inside of the apparent horizon. 
Therefore, if an apparent horizon appears, we may excise the interior of the apparent horizon to 
get rid of the physical singularities from the calculation domain. This idea (black hole excision) 
was proposed by W. G. Unruh [70], and had become a common practice to deal with physical 
singularities in numerical relativity. 

To perform black hole excision, one needs to locate the apparent horizon. 

5.5.1 Smoothness of Apparent Horizons in the Braneworld 

In the braneworld, the Israel’s junction condition at the brane essentially imposes cusp con¬ 
ditions to certain metric functions. This raises the question of whether apparent horizons will 
be non-smooth across the brane. By the discussion of the smoothness of apparent horizons in 
Sec. 2.2.2, the apparent horizons in the braneworld is smooth under the perpendicular gauge 
(2.7). Under gauge g rz | = 0 (perpendicular gauge), the smoothness can be simply expressed as 

(dr/dz) | = 0, where the derivative is evaluated along the apparent horizons. 

5.5.2 Apparent Horizon Finder 

In this subsection we introduce a method to locate apparent horizons [64]. 

We introduce polar coordinates (p, x) via r = p sin y r z = 1 + P cos x (length dimensions are in 
the unit of £). In this coordinate system, the symmetric axis r = 0 is y = 0, and the brane z = 1 
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is y = 7r/2. On a t = constant hypersurface, we define function 

/ = p- e(x), 


(5.29) 


and let the apparent horizon be the one with / = 0. The form of g is going to be determined by the 
apparent horizon equation (2.14). The spacelike outpointing unit vector normal to / = constant 
surfaces is now 


s 


Ot 


P a 

yJpPpp' 


where p a 


ladpf. 


(5.30) 


Substituting s“ into (2.14), we get a second order ordinary differential equation (ODE) of g with 
respect to y. The boundary condition at y = 0 (the symmetric axis r = 0 ) is dp/dy = 0 ; the 
boundary condition at y = 7r/2 (the brane) is dp/dy = 0 under the perpendicular gauge at the 
brane (fj rz = f] tz = 0 ), because of the Z 2 symmetry of the braneworld and the smoothness of the 
apparent horizon. Apparent horizon is obtained via solving this second order ODE subject to the 
boundary conditions. Numerically, we use a version of the shooting method to obtain the apparent 
horizon. There are multiple ways to implement shooting methods. The basic idea is to start the 
trajectory (implied by the ODE) at certain initial point (y = 0), and then solving the ODE subject 
to the initial conditions yields a value of dp/dy at the final point (y = 7 t/ 2 ) of the trajectory. We 
then adjust the inital point accordingly until dp/dy = 0 at final point. One version is implemented 
as the following. At y = 0, we pick up certain value of z( 0) p (which is g {y = 0)) as the initial value 
of z and solve the ODE to obtain the trajectory to y = tt/2. Not losing generality, let us assume 
the value of dp/dy at y = 7 t /2 is positive. Then we try different initial value of z(0) until we find 
an initial value z( 0) n such that dp/dy is negative at y = tt/2. Now we have found a bracket of the 
initial guesses: (z(0) p , t(0) ra ). Then we can use the following binary search to find the apparent 
horizon. We use z_p for z(0) p , z_n for z(0) n and Rp for dp/dy in the pseudo code. 


eps = pre-set small value 
Rp = 10 * eps 
do while ( abs(Rp) > eps ) 
z_m = (z_p + z_n) / 2 

solve ODE to get Rp at chi = pi/2, by initial value z = z_m 
if (Rp > 0) then 
z_p = z_m 

else if (Rp < 0) then 
z_n = z_m 
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end if 
end do 
z(0) = z_m 

In performing shooting method above, we need to solve the ODE, for which we used the Runge- 
Kntta method of second order accuracy. 

5.5.3 Dissipation at the Excision Boundary 

To perform the excision, which is to ignore the interior of the apparent horizon during numerical 
evolutions, is conceptually simple and neat. However, technically it is tricky and probably messy 
to deal with the excision boundary. To ignore the interior of the apparent horizon, we only evolve 
the exterior, and the boundary condition at the excision boundary is “no boundary condition”. 
Instabilities and noises tend to happen at the excision boundary, which should be removed by 
certain numerical dissipation. One way to do it, is to reconstruct the interior from the exterior, 
such that the reconstructed functions are smooth across the excision boundary. Then we apply 
Kreiss-Oliger dissipation to remove the noises. This process, while sounding trivial, is very difficult 
to realize. There is no good way to reconstruct smooth functions across the excision boundary. 
This method is able to evolve the spacetime with excision for “a while”, so its success is limited. 
Another way, which is the best among the methods we tried, is to use the one-side dissipation at the 
excision boundary, which was developed by Calabrese et al. [72], and adopted in PAMR/AMRD 
[65]. This method significantly improves the performance. 

5.6 Tests and the Validation of the Numerical Scheme 

As pointed out in Sec. 1.6.2, the tests after the numerical results are obtained, are essential to 
make sure the results are actually numerical solutions rather than numerical artifacts. To this end, 
the independent residual tests need to be performed. For systems with constraints, convergence 
tests for constraint residuals need to be performed as well. 

For GR, one can either perform independent residual tests together with constraint convergence 
tests, or perform convergence tests for residuals obtained from a different formalism of GR. Given 
generalized harmonic formalism (5.5) was used to obtain the solutions, where the conformally 
transformed metric and source functions were fundamental variables, now we perform the conver¬ 
gence tests of the residuals obtained from the original Einstein’s equations (5.1). i.e. the numerical 
solutions are obtained via the generalized harmonic formalism with the conformally transformed 
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metric as the fundamental variables. Then the solutions are transformed back to obtain the phys¬ 
ical metric, and are then substituted into the original Einstein’s equations to get residuals. These 
residuals should converge to zero at the expected order (which is of order 2 in our simulations), if 
the numerical solutions are obtained. The tests are shown in Fig. 5.3. 

Beyond this, we also performed independent residuals, which are perfect thus omitted here. 
Beyond these, we also performed convergence tests for the residuals of the constraints, which are 
shown in Fig. 5.3. 

5.7 The Numerical Solution 

The scheme can be used to study a wide range of dynamical processes, such as critical phe¬ 
nomena, the evolution problem in cosmology, gravitational wave from collapse, etc. But we will 
limit our attention to the end states of gravitational collapse at this time. 

We performed a series of simulations from the initial data metric (4.26) with the initial matter 
field as 

$ = g/ • exp [-(r - cc 0 ) 2 /(t^] , (5.31) 

from different initial data families. Within each family, only amplitude srf changes, while Gaussian 
parameters (ay and xq ), ay (the parameter in the “direct specification” function / in (4.11), which 
is used in eq. (4.31)) and the compactification parameters (ro and zo ), are fixed. 

5.7.1 The Evolution Process and Apparently Stationary State 

The initial data represents a localized Gaussian pulse. Since the initial data is time symmetric, 
the pulse evolves into two pulses: the ingoing pulse and the outgoing pulse. For weak data, the 
ingoing pulse is bounced back from r = 0 to travel outwards, which is the same phenomena as that 
in GR. The unique phenomena in RSII, is the interaction between the brane and the bulk, which 
mainly appears as the energy leaking into the bulk from the brane. Please refer to Fig. 5.4 and 
5.5. 

Sufficiently strong data will lead to BHs. The spacetime with BH can continue to evolve using 
BH excision techniques. The properties of the BHs are studied via apparent horizons in the bulk. 
Apparent horizon is generally different from event horizon. However, at the end of an evolution, 
the system reaches its stationary state and its apparent horizon coincides with its event horizon 
[89]. In the braneworld, when matter is absent around the horizon, the intersection of the bulk 


125 




5. 7. The Numerical Solution 


event horizon with the brane, is actually a well-defined event horizon on the brane as proved in 
Sec. 2.3, therefore observable, in principle. 

We monitor r a , Abuik and C 5 during an evolution, where r a is the areal radius of the intersection 
of the apparent horizon with the brane, Abuik is the area of the apparent horizon in the bulk, and 
C 5 is the length of the circumstance of the horizon (restricted on the r — z plane) in the bulk. For 
the simulations that can reach their end states or their apparently stationary state (to be defined 
below), the quantities reach the values that are almost constants. Please refer to Fig. 5.6 and 
Fig. 5.7 for an evolution which produces a BH of medium size, Fig. 5.8 for an evolution which 
produces a BH at a smaller size, Fig. 5.9 for an evolution which produces a small BH and Fig. 5.10 
for an evolution which produces a large BH. Fig. 5.4 and 5.5 show the evolutions from two initial 
data families, monitored by the Kretschmann scalar R^uapR^ 0 ^■ 

To obtain the end state of a system, it is necessary to let the evolution continue sufficiently 
long so that the system settles down. Even so, it is non-trivial to recognize the end state in a given 
simulation, due to coordinate effects. When a system reaches its end state, the system has a Killing 
vector that is asymptotically timelike, which corresponds to the time translational symmetry. If 
the slicing condition is such that the t = constant slices are Lie-dragged by the Killing vector, it is 
easy to overcome the spatial coordinate effects, by embedding the apparent horizons into a fixed 
space, such as the background space. During a particular evolution, the apparent horizon evolves 
with time, but its embedding in the fixed space will evolve into a time-independent state, if the 
slicing condition is adapted according to the Killing vector. 

However, it is non-trivial to impose such slicing conditions. Such a slicing condition is currently 
not imposed in our numerical schemes 16 . As a result, we do not have a time-independent state at 
the end of the evolutions. In the simulations, instead, we can obtain apparently stationary states — 
the states with BHs whose horizons (embedded into a fixed space) appear to be stationary for a 
finite time. By this definition, the apparently stationary state that stays stationary for infinite 
time is actually the end state. 

To recognize the apparently stationary state of a system, we monitor the evolution of the 
apparent horizon by its embedding into the “vacuum” background 

ds 2 = ^ — df 2 T df 2 + f 2 (d # 2 + sin 2 ddq ) 2 ) + d z 2 ^. (5.32) 

16 In the literature, there are the so-called symmetry seeking coordinate conditions [ 118 - 120 ] towards the conditions 
which can evolve into the coordinate configurations such that the time coordinate t is adapted to the Killing vector 
associated with the time translational symmetry of a stationary system. It will be our future work to implement 
such slicing conditions in our code. 
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Figure 5.4: The snapshots of four evolutions resulting in BHs with different sizes, monitored via 
Kretschmann scalar R^uafiR^' /a P. Each small panel represents an instant of the space, where the 
horizontal and vertical axes are r and z axes expressed in compactified coordinates r/(r + ro) and 
(z — £)/(z — i + zo), therefore the bottom of each panel is z = £, the brane. The complete space 
in r direction is shown. In z direction, only the part with interesting dynamics is shown. The 
evolutions are from the family with xg = 2,oy = 0.2, cr, = 0.4, and rg = zo = 2 are chosen. They 
produce no-BH, BH with r a = 0.29£, 0.61T, 3.78£ from stf = 0.04,0.15,0.24,0.49 respectively. The 
two smaller BHs are apparently stationary states, and the largest BH is not since the configuration 
has not settled down to stationary state after evolving for a long time, and eventually the code 
crashes. In general it is harder for an evolution to settle down if the excision surface is going across 
the “wiggling” regions. This appears to be a technical issue related to the black hole excision at 
relatively coarse resolutions. The black ellipses in the figures are the excision surfaces inside the 
apparent horizons. The evolutions clearly show how the energy flows from the brane into the bulk, 
and flow from the exterior towards the symmetric axis. Part of the energy is captured by black 
holes, the remaining energy continue to flow towards the “far end” of the bulk. 
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Figure 5.5: The snapshots of three evolutions resulting in BHs with different sizes, monitored 
via Kretschmann scalar R^vapR^ 1 ' 01 ^. The evolutions are from the family with xq = l,oy = 
0.1, cr 2 = 0.2, and ro = Zq = 1 are chosen. They produce no-BH, BH with r a = 0.22£, 2.9^ from 
£/ = 0.03,0.16,0.54 respectively. The smaller BH is a apparently stationary state, and the larger 
BH is not since the configuration is not settled down to stationary state after evolving for a long 
time. Again, presumbly this is a technical issue related to the black hole excision. 
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Here a bar (“) is used to emphasize that this spacetime is fixed, and the z > 1 portion is actually 
the background of the braneworld spacetime. The embedding is demonstrated by Fig. 2.3 (also 
eq. (2.40) and eq. (2.41)). The apparently stationary state appears as an “accumulating” curve in 
the embedding plot. Please refer to Fig. 5.6, Fig. 5.7, Fig. 5.8 and Fig. 5.9 as examples where the 
apparently stationary states are obtained. The processes of the settlement into time-independent 
states, show that the portion of the apparent horizons that is close to the brane, gets settled first, 
when the remaining portion might be still dynamical. 

Similar to apparent horizon, the existance of apparently stationary state is not guaranteed, and 
its relation with end state is not clear. In our simulations, as it turns out, apparently stationary 
states can be easily obtained by long term evolutions, as long as the apparent horizons do not cross 
the regions with interesting dynamics (the “wiggling” regions), which can be realized by properly 
choosing the initial data such that the wiggles either finally travel away from the apparent horizons, 
or are captured by the black holes. Furthermore, the plots of the quantities of apparently stationary 
states (such as the plots of -4b u ik-versus-r a and C 5 -versus-r a shown by Fig. 5.11), generated from 
the evolutions of the initial data profiles from different families , exhibite certain trends, whilst 
the same plots with horizons that are not apparently stationary state do not have trends. Also, 
the Abuik-versus-r a plot (the upper panel of Fig. 5.11) agrees perfectly with that obtained in the 
static system studied by Figueras-Wiseman in [20]. Therefore, we conjecture that the apparently 
stationary state is close to the end state, and we use apparently stationary state to approximate 
the end state. 

5.7.2 Black Holes as the Result of Gravitational Collapse 

For the BHs as apparently stationary states of gravitational collapse, we focus on the following 
aspects: the topology, the size and the shape. 

For the topology, one can see that the BHs appear to be localized on the brane with finite 
extension into the bulk. 

For the sizes, please refer to the results of all the simulations we performed, which are shown by 
Fig. 5.11, and table 5.1 for the results from some selected simulations. We did not try extremely 
hard to find the largest/smallest BHs possible. But within the simulations we performed, we 
obtained BHs with 

r a £ (0.04^,19.6^). (5-33) 

At the end of an evolution, the matter has either fallen into the BH, or escaped to infinity, 
which makes the brane tension be the only content associated with the brane. The strength of the 
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AH by Coordinates AH by Embedding 




(a) AH evolution by coordinates 


(b) AH evolution by embedding 



(c) The evolution of massless scalar field f l> 


Figure 5.6: This figure and Fig. 5.7 show an evolution that produces apparently stationary state as 
medium BH with (r a , _4buik, C 5 ) = (0.601,2.785,3.164), using parameters (£/,Xo,cr r ,a z ,ro,Zo) = 
(0.24,2,0.2,0.4,2,2) which defines the initial data and the compactification. (a) is the evolution 
of the apparent horizon, where r and z are coordinates. The black dashed line is the first apparent 
horizon that appears during the evolution, and the red line is the apparently stationary state 
of this specific evolution of gravitational collapse. Other lines are apparent horizons at intermediate 
instants, which color changes continuously from black to red. To better study the evolution without 
coordinate distortion effects, we embed each apparent horizon into the background spacetime: 
ds 2 = |j- ^ — df 2 + df 2 + r 2 (d 0 2 + sin 2 9d <^ 2 ) + d z 2 ^j . The embedding is demonstrated by Fig. 2.3. 

The embedding plot is shown in (b). This graph shows more clearly that the apparently stationary 
state is obtained: the shape of the apparent horizon has settled down, since the coordinate time 
t ~ 5.5. (c) shows the evolution of the massless scalar field that only lives on the brane. The initial 
profile is a Gaussian pulse, which splits into ingoing and outgoing branches. The ingoing branch 
is gradually captured by the BH. 
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0.62 



(c) The evolution of O5 (d) The minimum of d 

Figure 5.7: This is the continuation of Fig. 5.6. Here we show the evolution of r a , .Abuik, C 5 and 
the minimum of lapse function d over the whole calculation domain. One can see the apparently 
stationary state is obtained by sub-figure (b) in Fig. 5.6, since t ~ 5.5. The quantities approaches 
quasi-constant since then. However, the quantities do not stay strictly at constants. Combining 
the plot of d, hnally the values at t = 5.8 are recorded as the data for the apparently stationary 
state, min(d) ~ 0.35 means the lapsing rate of coordinate time is comparable with that of proper 
time. 
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AH by Coordinates AH by Embedding 




(a) AH evolution by coordinates (b) AH evolution by embedding 


0.24 

1 1 1 1 1 

1 1 1 1 

0.2 

111111111 






0.2 



0.15 


0.18 





0.16 



| 0.1 


0.14 



0.05 


0.12 

0.1 



0 



(c) The evolution of r a 


(d) The evolution of Abulk 
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(e) The evolution of C5 


(f) The minimum of d 


Figure 5.8: This figure shows an evolution that produces apparently stationary state as 
small BH with (r a , *4buik, C 5 ) = (0.221,0.184,1.311), using parameters (<£/, xo, oy> <Xz,ro> ^o) — 
(0.16,1,0.1,0.2,1,1). 
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AH by Coordinates 


AH by Embedding 




(a) AH evolution by coordinates 


(b) AH evolution by embedding 
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Figure 5.9: This figure shows an evolution that produces apparently stationary state as small 
BH with (r a , ylbuik, C 5 ) = (0.0432,0.00155,0.269), using parameters (^/, .To, ay, cr z , ?'o, zo) = 
(0.08,0.5,0.1,0.2,0.5,0.5). 
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AH by Coordinates 
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(a) AH evolution by coordinates 


AH by Embedding 
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(b) AH evolution by embedding 
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Figure 5.10: This figure shows an evolution that produces large BH as apparently station¬ 
ary state with (r a . Abulk, C$) = (15.1,2866,14.9), using parameters (jzf, xq, <x r , oy, ro, Zq) = 
(1.05,2,0.5,1,2,2). This is the case that data is so strong that all the interesting dynamics is 
captured by the BH. Note the plots of r a and .Abulk are noisy, which we can not explain. This 
might be caused by the KO dissipation—the dissipation appears to have stronger effects at larger 
spatial coordinates, and its perturbing effect on larger horizon is also stronger. As a result, the 
energy of the BH gradually decreases, which causes the decrease in Abulk- 
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Table 5.1: BHs produced from different initial data from different families. The initial data profile 
for the massless scalar field on the brane is $ = srf ■ exp [—(r — xo ) 2 /<t 2 ] • °z is the parameter to 
set up metric functions via “direct specification” eq. (4.11), which is then substituted into (4.31). 
The spatial metric for initial data is eq. (4.26). r o and zo are simply compactification parameters 
defined in eq. (5.9). 


brane tension is proportional to l/£, therefore invisible to small BHs whose size ?’ a -C £■ These 
BHs will be asymptotically 5D Schwarzschild, therefore Abuik = 27 t 2 r a and C 5 = 27 r?’ a . Please 
refer to Fig. 5.11. 

For the shape of large BHs, please refer to Sec. 5.7.3. 

5.7.3 The Relation with Black Strings 

The black string solution is ds 2 = (h a bdx a dx b + dz 2 ), where a = 0,1,2,3 and x a stands 
for a coordinate other than z, and h a t, is a BH solution of vacuum Einstein’s equations in 4D, 
which does not depend on coordinate z. This could be called black cone instead of black string, 
if we had considered the intrinsic geometry of the horizon of black strings. It is named string , 
in the sense that the [h a bdx a dx h + dz 2 ) part is a string. The “string” shape can be revealed by 
embedding the horizon into the background spacetime (5.32). Therefore, the embedding of the 
BHs into this background, can give a direct comparison with the black strings. The embedding of 
the BHs (represented by the apparently stationary states) of different sizes is shown in Fig. 5.12. 
Please refer to the caption of the figure for more details. 

Fig. 5.12 shows that small BHs are almost spherical. The BHs gradually change into a cigar- 
shape as the size increases, which suggests to call those BHs, black cigars , as first suggested in [ 6 ]. 
As the size of BH increases, the portion of the horizon that is close to the brane gradually changes 
into a black string. Also, the f 2 /z 2 factor in the background metric (5.32), contributes to Abuik as 
£ 3 /z 3 , which means the contribution from large z region become relatively negligible. As a result, 
for large BHs, only the portion that is close to the brane actually contributes to Abuik, regardless 
the behaviours in the large z region. This explains why the relation Abuik-versus-r a for large BHs 
(shown in the upper panel of Fig. 5.11), behaves as that of black strings. On the other hand, this 
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Figure 5.11: The shape of BHs with all sizes. Small BHs are asymptotically 5D Schwarzschild. The 
area-radius relation for large BHs is that of black string, which is also 4D Schwarzschild according 
to AdSs/CFT 4 [20, 22]. The area-radius relation is consistent with the one obtained from a static 
problem by Figueras-Wiseman. C$ for large BHs follows C 5 = 41og(£ • r a ) [23], where C = 2.71 is 
the best fit of our data. 
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r/{c + r a) 



Figure 5.12: The shape of BHs with all sizes shown by embedding into the background spacetime 
(5.32). The horizons are scaled by a factor of(c + r a ). Note both coordinates f and (z — 1) are 
scaled by the same factor, therefore the shape is not affected by the scaling. The scaling is to 
bring BHs with very different sizes (here r a ~ (0.04£, 15 £)) onto comparable plotting scale. Here 
the scaling parameter on the left panel is c = 0.3, and the right panel corresponds to c = 0. The 
left panel still carries the information of sizes, which emphasizes on the uniqueness feature. The 
right panel emphasizes on how the shape changes with size, which can also be read-off from the 
left panel. 

Small BHs are asymptotically 5D Schwarzschild, which are almost spherical. They then change 
into cigars (seen from the embedding point of view) as the size increases, which suggest to call 
them black cigars as firstly suggested in [6]. 

More importantly, the BHs are apparently stationary states from different families. For family 
(i), the parameters specifying the family are (xo, cr r , cr z , rg, zo) = (0.5, 0.1, 0.2, 0.5, 0.5). These 
parameters for other families are: (ii) (0,0.3,0.3,1,1); (iii) (1, 0.1,0.2,1,1); (iv) (2,0.2,0.4, 2,2); 
(v) (2,0.5,1, 2,2). From the left panel we see that BHs obtained from different families with the 
same size, almost agree with one another. This indicates the detail of initial data is lost, and the 
solution is unique. The right panel shows the relative extension into the bulk increases with the 
size, and the portion (of the horizons) that is close to the brane looks more and more like black 
string as the size increases. 
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also means that the relation of *4buik-versus-r a can not reveal the difference between the BHs and 
black strings. 

Therefore, we study C 5 (the length of the circumstance of the horizon restricted on the r — z 
plane) versus r a , since C 5 is infinite for black strings, but finite for BHs. Furthermore, Fig. 5.12 
shows that the relative extension into the bulk (as seen from the point of view of the background 
metric (5.32)), increases with the size on the brane. It is not very clear whether the relative 
extension has an upper limit. On the other hand, by the uniqueness of the BH solutions that is 
going to be studied in Sec. 5.7.4, and the comparison with Figueras-Wiseman solution that is going 
to be shown in Sec. 5.7.5, there is strong indication that our solutions generated from evolution 
systems, are the same as those were obtained by Figueras-Wiseman in their static systems. For 
the Figueras-Wiseman solution, as shown in the upper figure of Fig. 5.14, the large BHs has a 
limiting shape which is the AdSs/CFT 4 solution [20, 21 ]. Generally, for configurations with fixed 
shape in the space (5.32), there is the following property for large r a 

C 5 = 41og(C • r a ), (5.34) 

regardless what the shape is, as long as the shape is fixed. This statement can be justified by direct 
numerical experiments for a few shapes. £ is determined by the specific shape. C 5 versus r a is 
plotted as the lower panel of Fig. 5.11, which supports the assumption that large BHs have the same 
shape , and the best fit of £ is 2.71. C ~ 2.8 was first independently found by Figueras-Wiseman 
[23], and eq. (5.34) was proposed by Toby Wiseman in [23]. 

Large BHs have the same shape (black cigar) in the background space (5.32), and black strings 
also appear as strings (rather than cones) in this space, whilst C^/r^ merely appears as log(r a )/r a 
(black pancake) which is not as great as a fixed shape (personal tastes), therefore it makes more 
sense to name BHs as black cigars as first suggested in [ 6 ], rather than black pancakes. 

5.7.4 The No-hair Feature 

We purposely performed the simulations from distinct initial data families. Because the result 
of a well-posed numerical simulation depends smoothly on its initial data, if only one family is 
considered (only one parameter srf in the initial data profile is changed), the quantities (such as 
-Abuik, C 5 and r a ) will depend smoothly on srf . Therefore, for a given family, relations between 
quantities such as Abuik-versus-r a , will emerge. For a different family, in principle the relation of 
A b uik-versus-r a might be different from the relation obtained from the previous family. Should 
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this happens, the BH solutions are not unique. In reality, however, the relations of .Abuik-versus-r a 
obtained from different families are the same: the .Abuik-versus-r a relation plotted using the data 
obtained from the evolutions from different initial data families is shown in the upper panel of 
Fig. 5.11, which appears to be one single curve. Similarly, C 5 -versus-r a is plotted as the lower 
panel of Fig. 5.11, which also appears to be one single curve, i.e. Fig. 5.11 shows that the BHs 
with the same sizes have the same areas and the same circumferences, regardless which families 
the results are generated from. This indicates that the shape of the horizon is solely determined by 
the size r a , regardless which initial data family the BHs are generated from. Therefore a no-hair 
theorem of the BH solution in RSII is suggested. In general, however, the BH solutions in AdS 
spacetimes may not be unique (see, e.g. [117]). Therefore, the uniqueness of the BHs is limited 
to the RSII spacetimes studied in this thesis—these spacetimes are axisymmetric without angular 
momentum and non-gravitational charges. In this situation, the shapes of the horizons are directly 
studied in Fig. 5.12, where one can see that the BHs with the same size produced from different 
families actually have almost the same shape, which is an indication that the detail of the initial 
data is lost in the final state, and the geometry of a BH is solely determined by its size. Therefore 
a no-hair theorem (the uniqueness of BH solution) is suggested to hold for BHs in RSII. 

5.7.5 The Comparison with Figueras-Wiseman Solution 

If the BH solutions in RSII are unique in axisymmetric spacetimes without angular momentum 
and non-gravitational charges, then the BH solutions should be the same, regardless how the BH 
solution is obtained. In particular, the BHs produced as the end states of evolutionary systems 
should agree with the ones obtained from the static problem studied in [20, 22]. Here we still use 
apparently stationary states to approximate end states. Following [20, 22], we plot our data of 
Abuik versus r a on top of the Figueras-Wiseman data in the upper panel of Fig. 5.11. The figure 
shows the agreement with Figueras-Wiseman solution as illustrated by Abuik versus r a is perfect. 

To better compare with their solution, we also embed the BHs into the space (5.32) with z < 1, 
which is what Figueras-Wiseman did in [20]. Following [20], the freedom of the embedding is fixed 
by mapping the r = 0 ends (i.e., the axis ends) of the horizons to the point (f, z) = (0,1) in (5.32), 
instead of mapping the z = 1 ends (i.e., the brane ends) to z = 1 in (5.32) as what we did in 
the above sections, i.e. instead of performing the embedding as Fig. 2.3, here we perform the 
embedding as Fig. 5.13. Please refer to Fig. 5.14 for the results. Fig. 5.14 shows the two results 
qualitatively agree. However, there are a few differences between these two figures: the largest BH 
meets the vertical axis at f « 0.468 for our data, but at f « 0.457 for Figueras-Wiseman data; there 
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(a) 


(b) 


Figure 5.13: The embedding of a closed surface. Fig. (a) is the closed surface in the physical 
spacetime, and Fig. (b) shows its embedding into the background spacetime (5.32) with z < 1. 
The freedom of the embedding is fixed by mapping D to D. 


exist line crossings (can be seen by zooming in) in our data but there is not in Figueras-Wiseman 
data. Together with the line crossings in Fig. 5.12, it implies that generally apparently stationary 
states are, close to but distinct from, end states. Here we emphasize that our solution for large 
BHs were obtained at a resolution that is effectively coarse at the horizon, thus the solution should 
be less reliable. Furthermore, as one can see from the processes of the settlements to apparently 
stationary states (Fig. 5.6, Fig. 5.7, Fig. 5.8 and Fig. 5.9), the portion of the apparent horizon 
that is close to the brane gets settled first, while the remaining portion (the portion that is far 
from the brane, and close to the symmetric axis) might be still dynamical. This makes the portion 
that is close to the brane more reliable than the portion that is far from the brane. Therefore, 
the embedding plot by fixing the brane end of the apparent horizon as done in Fig. 5.12, is more 
reliable than the embedding plot by fixing the axis end of the apparent horizon as done in Fig. 5.14. 

5.7.6 Brane Energy 

In this section we focus on the physics on the brane, by studying the quantities obtained from 
the reduced metric (h^) on the brane. The rationale is that the brane is all one can directly 
observe (while the bulk is invisible). We compare quantities obtained on the brane, with the same 
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(a) Figueras-Wiseman solution 



Figure 5.14: The comparison with the Figueras-Wiseman solution via embedding into (5.32), 
with z < 1. The freedom of the embedding is fixed by mapping the r = 0 ends of the horizons to 
z S3 1 in (5.32). Note, figure (a) is from [20], and the labels (r, z) should be understood as (r, z). 
The large BHs have a limiting shape which is the AdSs/CFT 4 solution [20, 21]. In (b), the largest 
BH is the one with r a = 15.1. 
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quantities obtained from 4D GR. We will study ADM mass and Hawking mass. As explained in 
Sec. 2.5, these masses are only defined at asymptotic spatial infinities. However, let us first study 
the ADM mass and Hawking mass in 4D GR, to obtain some insights. 

To compare with the braneworld where a spherical symmetry is present on the brane, we first 
consider gravitational collapse under spherically symmetric configuration in 4D GR. The initial 
data is time symmetric Gaussian (massless scalar field), centred at r = 1. The pulse splits into 
ingoing and outgoing branches. For data that is not strong enough to produce BH, the ingoing 
branch is bounced back and then travel outwards. For the ADM mass and Hawking mass in 4D 
GR, please refer to Fig. 5.15. The figure shows that: 

(1) the energy is conserved (remains to be constant at spatial infinity); 

(2) the energy is monotonic with radius r, which means the energy “density” (not well-defined in 
general) is positive; 

(3) there is a “stair” between two separate matter pulses—this is due to Birkhoff’s theorem: 
vacuum solutions in spherical symmetry are unique (which are Schwarzschild’s solutions). 
Therefore, in the vacuum gap between the two pulses, the solution is should be Schwarzschild, 
where the ADM/Hawking energy is well-defined as well. Thus Hawking mass is constant within 
this gap, while ADM mass is approximately constant. 

For the ADM mass and Hawking mass on the brane, please refer to Fig. 5.16, where one can 
see that the above mentioned features (2) and (3) are lost. Feature (1) still holds—this is because 
the masses are defined as the limit at the spatial infinity utilizing the time translational symmetry. 
Or in another word, this is because the dynamics happening locally can not propagate to the 
spatial infinity within finite time. However, it is expected that feature (1) does not to hold in the 
braneworld. This is because, there is energy exchange between the bulk and the brane, and the 
simulations show that the “wiggles” are moving from the brane into the bulk. Therefore we need 
a quantity that can describe this phenomena. 

It turns out the energy defined by equation (2.69) can indeed fulfill this task. Please refer to 
Fig. 5.17. This energy has feature (2) and feature (3) mentioned above. The energy at spatial 
infinity, on the other hand, changes with time. Fig. 5.17(b) shows the change with time: the 
brane loses energy in majority of time (which agrees with the phenomena we observed from the 
simulations we performed), but it gains some energy when the incoming pulse gets bounced back 
at the origin. 


142 




5. 7. The Numerical Solution 


GR ADM Energy GR Hawking Energy 




Figure 5.15: Energies in 4D GR. Each line is plotted at an instant (i.e., constant coordinate time). 
The black dashed line is at the earliest instant, and the red line is at the last instant. Other 
lines’ colour changes gradually from black to red, with respect to coordinate time. Each line only 
has one colour since it stands for one instant. 


Brane ADM “Energy” 



Brane Hawking “Energy” 



Figure 5.16: “Energies” on Brane. Each line is plotted at an instant (i.e., constant coordinate 
time). The black dashed line is at the earliest instant, and the red line is at the last instant. 
Other lines’ colour changes gradually from black to red, with respect to coordinate time. Each line 
only has one colour since it stands for one instant. 
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Figure 5.17: This figure shows the brane energy defined by equation (2.69). In (a), each line 
is plotted at an instant (i.e., constant coordinate time). The black dashed line is at the earliest 
instant, and the red line is at the last instant. Other lines’ colour changes gradually from 
black to red, with respect to coordinate time. Each line only has one colour since it stands for one 
instant, (b) shows the energy at spatial infinity changing with coordinate time. 
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Chapter 6 


Conclusion and Future Work 

The basic idea of braneworld scenarios is that our observable universe could be a 3+1 dimen¬ 
sional surface (the “brane”) embedded in a higher dimensional spacetime (the “bulk”). The single 
brane scenario constructed by Randall and Sundrum (also known as RS-braneworld II), is what 
we study in this thesis project. The basic setup is as follows: the single brane (the observable 
universe) is embedded in the bulk with one extra dimension of infinite size. The matter is strictly 
trapped on the brane while gravity is free to access the bulk. The bulk is therefore “empty”, but 
has an assumed negative cosmological constant. The bulk has a Z-± symmetry with respect to 
the brane, and, the brane has a tension which enable fine-tuning to any equivalent cosmological 
constant on the brane. General relativity (GR) is recovered on brane at low energies, but the 
brane dynamics can be rather different from GR at high energies. This latter regime is the focus 
of our research. Specifically, we study the dynamical process of black hole formation as a result of 
gravitational collapse of massless scalar field. 

6.1 Conclusion 

We have achieved the following: in terms of developing the machinery, we discover/develop/in¬ 
vent several novel facts, formalisms and techniques regarding NR and braneworld. The regularity 
problem in previous NR simulations in axisymmetric (and spherically symmetric) spacetime, is 
actually associated with neither coordinate systems nor the machine precision. The numerical cal¬ 
culation is regular in any coordinates, provided the fundamental variables (used in simulation) are 
regular, and their asymptotic behaviour at the vicinity of the axis (or origin) is compatible with the 
finite difference scheme. Generalized harmonic (GH) and BSSN formalisms for general relativity 
are developed to make them suitable for simulations in non-Cartesian coordinate under non-fiat 
background. A conformal function of the metric is included into the formalism to simulate the 
braneworld. The usual constraint damping term used in GH, can not control the severe constraint 
violation in braneworld. The violation is cured by imposing the constraints properly. In solving 
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elliptic equations (Hamiltonian constraint, for instance), using functions to carry the asymptotic 
behaviour at spatial infinities could be crucial. The delta-function matter can be simulated by 
“integrating out” the delta function and the brane content can be encoded in Israel’s junction 
condition. 

On the physics side, we perform the first numerical study of gravitational collapse in braneworld 
within the framework of the single brane model by Randall-Sundrum (RSII). The scheme is capable 
of obtaining apparently stationary states as the results of gravitational collapse, which are BHs 
localized on the brane, with finite extension into the bulk. The extension changes from sphere 
to flattened pancake (or from sphere to cigar, from the embedding point of view) as the size 
of BH increases. There is strong evidence that the detail of initial data is lost in the resulting 
BH, therefore no-hair theorem of BH (uniqueness of BH solution) is suggested to hold in RSII 
spacetimes that are studied in this thesis—these spacetimes are axisymmetric without angular 
momentum and lion-gravitational charges. In particular, the BHs we obtained as the apparently 
stationary states from the dynamical system, are consistent with the ones previously obtained from 
a static problem by Figueras and Wiseman. We also obtained some results in closed form without 
numerical computation: the smoothness of apparent horizon across the brane under perpendicular 
gauge, the well-definedness of event horizon on the brane, and the equality of ADM mass of the 
brane with the total mass of braneworld. 


6.2 Future Work 

There are many potential research directions/projects suggested by this project, some of which 
will become our future work. One such direction is to derive the mass-area relation of the (apparent) 
horizons in closed form. Another project is to look at the examples related to holographic principle: 
the equality of bulk energy and the ADM/Hawking energy calculated from the brane geometry 
which we have proved to hold for two classes of spacetimes. We would like to examine whether 
the relation holds in general space, and study the asymptotical configurations at spatial infinities. 

As for this work per se, it can be expended/upgraded as follows 

(1) Apparently stationary states were obtained, but it is not clear whether end states were ob¬ 
tained. The method to identify end state needs to be developed. Some attempts would be to 
implement the slicing conditions in [118-120]. 

(2) The current slicing condition and coordinate gauges can perform a wide range of initial data 
and yield apparently stationary states. However, for certain initial data, the gauges we used 
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(Harmonic gauge, lapse driven gauge, damping wave gauge) can not yield a simulation that 
last “forever”. Therefore we need to improve the performance by studying coordinate gauges. 

(3) There is slight inconsistency in the initial stages of the evolution (which fade out during 
evolution), which might be due to the truncation error, or inconsistency in initial data. Further 
study is needed. 

(4) The energy of the brane is needed to describe the interaction between the brane and the 
bulk. The brane energy we obtained is not conceptually satisfying, although the qualitative 
behaviour is surprisingly good. The energy aspect needs more study. 

(5) The current code can run in a parallel environment since it is based on AMRD/PAMR. How¬ 
ever, we used a single CPU to do all the simulations since shooting method were employed to 
search for the apparent horizon. One improvement is to use flow method to search for apparent 
horizon, which can be used in a parallel environment. 

For the same reason, although the code is capable of using the adaptive mesh refinement 
(AMR) that is built-in in AMRD/PAMR, because of the shooting method used in finding 
apparent horizon, we performed the simulations under unigrid. Hopefully we can turn on the 
AMR feature when the shooting method is replaced by flow method. In this way we can also 
improve the accuracy of the result. Here we emphasize that resolution is important since our 
calculation of large BHs resulted in apparent horizon at z ~ 30f, where the effective resolution 
was quite low and should be improved. 

Also, we could try another way (which I called BSSN-like method) to implement the generalized 
harmonic formalism as discussed in section 5.3. 

(6) Generalized harmonic formalism of GR was used in this work. We could also try other formal¬ 
ism such as BSSN, where it is more natural to impose gauges. 
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Appendix A 


Generalized BSSN Formalism 


The original BSSN [45, 46] formalism is only applicable to 4D spacetime with Cartesian co¬ 
ordinate. It was generalized into non-Cartesian coordinate in [36, 97] for 4D spacetimes with 
flat backgrounds. We will generalize BSSN into arbitrary dimension, general background, general 
coordinate and general conformal function (without requiring 7 = 7 , see below), yet free from 
irregularity issue. This work is essential to make BSSN applicable to higher dimensional space- 
time with non-trivial background, such as braneworld simulation. The derivation in [36] is largely 
borrowed in this section. 

Fundamental variables to evolve: ip, K, "fij , A^, T®, which are defined as 


7 ij = e w 7 ij, 

(A.l) 

A t j — ^ Ky d _ ^7) 

(A.2) 

Aij = e qv Aij, 

(A.3) 

P _ ~j k Ad-i)f i. k _ («*-i)r< \ , 

(A.4) 


where 7 , the determinant of 7 ij, is set to a time-independent function, such as 7 (the determinant 
of the background spatial metric). The bar-e d quantities are associated with the time-independent 
background metric. 

Let’s derive the evolution equations for these quantities. 

The starting points are the evolution equations for I\ a p and ^ a p in ADM-York formalism: 
eq. (1.15) and (1.16). They are written below as 


Cm — eD a Dpa -f- ot 


-e < rf : + KK a p 


2 Ka/xK^p + £ kd 


7'>/3 


S_+ep\ 

d -2 ) 


(A.5) 


f^m7 a/3 — 2o:A a( g, 


(A. 6 ) 
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from which we can derive 

£ m ln 7 1/2 = ^7 y £ m 7 ij = ~aK. (A.7) 

The equation of motion for K is 

C m K = + K u C m -> iJ 

= eD'DiO + a [i„I« + ^* k„ ((£§„ - , (A.8) 

where the Hamiltonian constraint is used. The equations of motion for ip, jij are derived to be 
[An ln 7 1 / 2 -£m ln 7 1/2 ] = 9 ( d _i) (-a^-^mln7 1/2 ) . (A.9) 

£m7t? = e~ w C m Jij - qe~ w ^ijC m p = -2aAy + ^ 7ij£ m lny 1 / 2 . (A.10) 

For T\ we have 

f ,: = 7 jfc ( (d - 1) f i ifc - (rf = -Dir + f j dj [ln(7/7) 1/2 ] ■ (A.ll) 

=7 d t r = -D j (dty») + dt^dj [ln(7/7) 1/2 ] • (A.12) 

where we have applied = i Q k In 7 (and ( d_ 1 )r J , k = 1 dk In 7 ), and dtDj — Djdt = 0 which 

is because < 9 t ( d-1 )p , k = 0. dt^ can be easily evaluated from (A. 10) via <9t7 y = —7* fe 7 ji 9t7fc;. 
The evolution of Ay needs the result of the evolution of Ay, which is 

I'm -1 ij —Tin 1 \ y — - (A £ m 7y T 7y T m /\ ) 

= (£ m A'y) TF + - ^^nK lm K lm , (A. 13) 

where TF means trace free. Therefore 

Ay = e~ qv ’C m A ij - qe~ qv A i:j C m ip 

■ / o , \ _ TF 

= e~ w eDiDja - ae (d_1) i?y + eak d ( Sy - 7y-^—y J 
+ aKAij - 2a An A 1 j + ^ ^ ^ A VJ L m lnqA/ 2 . (A.14) 

( d_ ^Rij needs to be rewritten in terms of the F* defined in (A.ll). By the same procedure for the 
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conformal transformation carried out in [32, 36], we have 

= ^R tj - q(d ~ 3 ) DiD j¥ > - !%■ D k D k <p + r{d ~ 3) (p i{p D jlp - % D k <p D k <p ). 

(A.15) 

Therefore 

= e - qv (v-^R - q(d - 2 )b k D k ip - q2 ^ d ~ ~ 2) D k ipD k ip\ . (A.16) 

Expressing ( d ~ 1 and in terms of E 1 , has been done in section 3.3 (or [32, 36]) 

(d “ 1) Rij = {d ~ 1] Rij ~ \i kl R k Dan - - D (i Yj) + f k C k tj - C k u C l kj . (A.17) 

Here we have defined C l - fc = ( d_ 1 )f* - k — ( d_ 1 )T* . fc , and f; = t^T-L Contracting (A.17) with y*- 7 , 
we can obtain the expression of R in terms of T*. 

In (A.14), DjDja needs to be expressed by its counter parts 

DiDja = D;Dja = DiDja — C k i jD k a 

= DiDja - | (DiaDjip + DjaD^p - y D k ipD k a^J , (A.18) 

where = ( d ~ 1 \' k i . — and = § ^<5 k iDj^> + 5 k jDpp — jijD k (p'j are used, and the 

latter can be obtained by repeating the relevant derivations in [36]. 

To perform numerical relativity, all the equations above need to be rewritten as partial deriva¬ 
tives with respect to t, which are related to C m by £ m = Cg t — Cp. We must be careful with the 
Lie derivatives of tensor densities with respect to /3. An object X is a tensor density of weight w, 
if X = tensor x y' 0 ’/ 2 . Its Lie derivative is 

CpXmlCpX} +wXd i /3 i . (A.19) 

L J w—0 

Let’s now figure out the weight of fundamental variables. Because of (A.l), we have = 
( 7 / 7 ) qid ~ 1) , therefore the weight of e v is where w is the weight of 7 . Similarly, the weight 

of 7 ij and A.jj is the weight of upstairs y *- 7 and A '- 7 is 737 - The value of w, can actually 

be determined as follows. 7/7 is a scalar, because its value does not change under coordinate 
transformation. Therefore, 7 = (7 • tensor), which implies w = 2. Therefore, e v , 7 ij and Aij are 
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all tensors 1 '. Knowing w = 2, it is then easy to derive 

jCp In 7 1 / 2 = -^£/37 = • • ■ = fd l di lny 1 / 2 + d,/?*. (A.20) 

27 

The formulae can then be rewritten in terms of coordinate derivatives by opening Cp ln^ 1 / 2 . Let’s 
take (A.9) as an example. Using m = dt — f3 and dtj = 0, we have 

d ' v = w^T) h K +« + K 1 ^+< A21 > 

where we have applied Cg t T = dt.T in a coordinate system where t is a coordinate, for tensor T. 
The Hamiltonian constraint is now 

k d p * \ (-e«-VR + , (A.22) 

where is of course replaced by its expression in equation (A.16). ~ means the equation is 

a constraint relation. The momentum constraint reads 

ek d Si ~ D,K - (DjAF + qi - d ~ l) D jip A\ + 

= ~ (CC 1 ^%'C) • (A.23) 

Summary 

The fundamental variables are ip, K, 7 ,;.,. Aij, T®. The equations of motion are (A. 8 , A.9, A.10, 
A.12, A.14). Constraints are the Hamiltonian constraint (A.22), momentum constraints (A.23), 
and the definition equations (A.1-A.4), which read 

7 ~ the pre-set time-independent function such as 7 , (A.24) 

7 ij Aii - 0, (A.25) 

I-'" A r '.h.)- ( A - 26 ) 

17 Quite a few authors counted the weights incorrectly—especially under Cartesian coordinate with flat background 
where 7 was set to 1, where these authors incorrectly assumed w = 0. Fortunately, the values of the weights per se 
do not matter in the final expressions in terms of partial derivatives. It is the relative weight that matters, e.g. one 
can keep w general (without substituting w = 2), therefore the weight of is • One can obtain (A.21) 

correctly. 
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Axisymmetry 

For simulations in axisymmetry—take 4D as an example—under cylindrical coordinates (t, p , <f>, z), 
as shown in Chap. 3, the such defined F* reduces to the results obtained by our Cartesian com¬ 
ponents method if the background 7 ^ is flat, therefore regular. When the background is not flat, 
the behaviour needs to be analyzed in a case-by-case basis. The conformal metric and conformal 
traceless extrinsic curvature are 



where local flatness has been applied. All the fundamental variables depend on (t. p , z ) only. 
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Extrinsic Curvature as Geodesics 
Deviation: C > 1 Case 


In this section, we consider the case that a (d — C) dimensional object E embeds in the d 
dimensional space M. Therefore C is the codimension. 

Let ... be C continuous normal vector fields of E, with unit length ( 1 ^e, ... , ^ c \, 

and are set to be mutually orthogonal, where the index I = 1,2,...,C. Each = ±1, 

depending on the spacelike/timelike nature of the corresponding dimension. 

The procedure to obtain equation (2.23) can be straightforwardly generalized to C > 1 case. 
For example, generalize equation (2.17) to 


lafi — 9a.fi 


(1) n a (1 V — — (C ^e ^ c \ia- c \ip. 


(B.l) 


Eventually the deviation of two geodesics produced by T“ £ E is derived to be 

c c 

(/) e = T a T p ^2 {I \ (/) n M {I) I< a0 . (B.2) 

7=1 1=1 

i.e. we have defined C tensors 

{I) K a0 = 7 / 7 /W (< 7 W) , where 7=1,..., C. (B.3) 

On the other hand, NVP takes the “change rate” of the direction of unit normal vector along 
E to serve as extrinsic curvature. When C > 1, the direction of C dimensional orthogonal space 
is characterized by the wedge form 

(1 W...A (C V. (B.4) 
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Therefore, under NVP, the extrinsic curvature is 18 

K Sli ...„ = 7 ,“V a ( (1) n„ A ... A (c) n,) , (B.5) 

which is a tensor of rank C + 1. As the form implies, _ v does not lie within E. However, we 
can prove that NVP is still equivalent to GEP. 

First, we define the notation 

ds = 7 /*V Q . (B.6) 

Then it is easy to show that 

Ks»p...v = ds A (2) rig A ... A A (2) ng A ... A (c Vi„ 

+ (1) n M A ( ' dtf^np ) A ... A (C W + (1 Vi M A (2) n^ A ... A (d ( s) (c \i^ . (B.7) 

i.e. Leibniz rule. And of course the operation ds should be excluded from the wedge’s. We use a 
bracket on <5 as d(S) to express this fact explicitly. 

Before we continue, let us introduce some notations 

c 

A T liV = g^u - 7 m ^ = (/)e ^n^riv. (B.8) 

7=1 

W^ v wJV M1/ -We( J W= ]T (J) e (J V J W- (B.9) 

Noticing that = P) e = constant along E, we have 

(d« (J W) {I) n v = 0, (B.10) 

which implies 

d<5 (J W = 7,$“V a (/ W = (^75 a V a (/) n M ) 5/ 

= ( 7 /^^%) ( 7 / + A/) = ( 7 * a V a «n M ) ( 7 / + , (B.ll) 

18 Or define the extrinsic curvature via the derivative of the wedge of the D — C tangent vectors of E. 
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where we have applied equation (B.10) in deriving the last equal sign. The above equation is just 

ds (I) n u = (I) K 6 „ + ^ {J W IJ) B S , (B.12) 

J=i, J// 

where we have defined 

(U) Bg = ( (J W*d 4 <%) ■ (B.13) 

Overall, the above procedure is nothing but projecting d$( I \i v into the tangent space of E and the 
orthogonal space of E. 

Now we can use the fact 

... A A ... A A ... = 0, (B.14) 

and use equation (B.12) to rewrite equation (B.7) as 

Ksup... v = ^ A ... A ^n v 

+ (1) n M A (2) K m A ... A (c) n„ + • • • + (1) n M A A ... A {C) K {5> . (B.15) 

Again, when doing wedge, S is not affected. The terms above are linearly independent, therefore 
Kg/id.-.v is a tensor whose coefficients of linearly independent tensors are ^K a p. In another word, 
there is a one-to-one correspondence between Kg and a set of B)K a p. Therefore, NVP extrinsic 
curvature Kgn.p...v is equivalent to C quantities < ' r ' l K l , M , which characterize extrinsic curvature under 
GEP. i.e. GEP and NVP are equivalent for any codimension. 
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